http://orcid.org/0000-0002-1769-1368
Figures
Abstract
Transcriptome deep sequencing is a powerful tool for exploring the genetic architecture of complex traits. Gene expression patterns may explain a high degree of the observed phenotypic differences in histochemical and metabolic parameters related to meat quality among different muscles. In this study, we sequenced by RNA-Seq the whole transcriptome of nine lamb muscles: Semimembranosus (SM), Semitendinosus (ST), Cranial gluteobiceps, Gluteus medius (GM), Rectus femoris, Supraspinatus (SS), Longissimus lumborum (LL), Adductor and Psoas major. Significant gene expression differences were detected between almost all pairwise comparisons, being more pronounced between SS and ST, SM and LL, and ST and GM. These differences can be explained in terms of ATPase and glycolytic activities, muscle fiber typing and oxidative score, clustering muscles as fast glycolytic, intermediate or slow oxidative. ST showed up-regulation of gene pathways related to carbohydrate metabolism, energy generation and protein turnover as expected from a fast white muscle. SS showed myosin isoforms typical of slow muscles and high expression of genes related to calcium homeostasis and vascularization. SM, LL and GM showed in general intermediate gene expression patterns. Several novel transcripts were detected, mostly related to muscle contraction and structure, oxidative metabolism, lipid metabolism and protein phosphorylation. Expression profiles were consistent with previous histochemical and metabolic characterization of these muscles. Up-regulation of ion transport genes may account for significant differences in water holding capacity. High expression of genes related to cell adhesion, cytoskeleton organization, extracellular matrix components and protein phosphorylation may be related to meat yellowness and lower tenderness scores. Differential expression of genes related to glycolytic activity and lactic acid generation among fast, intermediate and slow muscles may explain the detected final meat pH differences. These results reveal new candidate genes associated with lamb meat quality, and give a deeper insight into the genetic architecture of these complex traits.
Citation: Armstrong E, Iriarte A, Nicolini P, De Los Santos J, Ithurralde J, Bielli A, et al. (2018) Comparison of transcriptomic landscapes of different lamb muscles using RNA-Seq. PLoS ONE 13(7): e0200732. https://doi.org/10.1371/journal.pone.0200732
Editor: Cristina Óvilo, INIA, SPAIN
Received: August 7, 2017; Accepted: July 2, 2018; Published: July 24, 2018
Copyright: © 2018 Armstrong et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Sequencing data can be accessed through NCBI GEO with accession number GSE112500.
Funding: Funding source: Comision Sectorial de Investigacion Cientifica (Universidad de la Republica, Uruguay). ID of the project: I+D 679 (year 2012). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
High and consistent quality are critical for the sheep meat industry, which could be achieved through improving breeding strategies and management procedures but also by identifying genes that can help to predict meat quality [1]. The latter is especially critical in sheep, as little is known about its genome and gene expression profiles compared to cattle, pig or poultry.
Meat quality involves several complex traits, difficult to improve and control. Transcriptome deep sequencing is a powerful tool for unraveling the genetic architecture of these traits. Gene expression patterns explain a high degree of the observed phenotypic variation, and thus the identification of genes and metabolic pathways involved could help to improve this product [2–4].
Previous studies of gene expression profiles in sheep muscle include comparisons between different breeds through RNA-Seq [5–7], and muscle microarray-based studies of callipyge vs. wild-type lambs [8, 9] and microRNAs [10]. Sun et al. (2016) found 960 genes differentially expressed (DE) in Longissimus dorsi muscle from two Chinese sheep breeds that differ greatly in meat production traits, mostly related to muscle growth, muscle cell differentiation, metabolic properties and cellular components [7]. In a similar study comparing muscle Biceps brachii in divergent sheep breeds, most DE genes were involved in structural and metabolic pathways [5, 6]. Differences in gene expression between callipyge and wild-type lambs related mostly to muscle fiber type conversion and muscle growth [8, 9]. Zhang et al. (2015) found 22 miRNA highly expressed in sheep Biceps femoris muscle, some of which were associated with meat production performance [10]. To our knowledge, there are no studies that compare general expression profiles of different sheep muscles under normal conditions, as has been done in other species [11].
In mouse, transcriptome deep sequencing of a fast muscle (Quadriceps) in comparison with a slow one (Soleus) showed pronounced differences in overall gene expression, especially in genes related to muscle contraction, carbohydrate metabolism and different types of non-coding RNA [11].
Mammalian skeletal muscles differ in their contractile and metabolic properties, mainly related to speed of contraction and resistance to fatigue. In a previous study, Ithurralde et al. (2015, 2017) studied fiber content and type among lamb muscles with different histochemical and meat quality characteristics [12, 13]. The authors demonstrated the existence of specific intermuscular differences in fiber typing and meat quality in sheep. How these properties relate to gene expression profiles and ultimately to meat quality is largely unknown, and this knowledge is of great interest for the meat industry. In this study we sequenced the whole transcriptome of nine different muscles that belong to highly valuable cuts in unrelated male lambs using RNA-Seq technology, in order to associate gene expression patterns with histochemical properties, metabolic parameters and meat quality characteristics. Our results indicate that transcription profiles differ between muscles, especially according to their metabolism and main fiber types, and that the detected gene expression differences may explain part of the muscle phenotypic variation observed.
Materials and methods
Sample collection, muscle fiber typing and meat quality measures
All animal experimental procedures were approved by the Animal Experimentation Committee (CHEA, Universidad de la República, Montevideo, Uruguay). Muscle samples were collected in a previous study by Ithurralde et al. (2015, 2017) [12, 13] from four unrelated Poll Dorset cross-bred lambs, raised on grasslands. Animals were acquired from different flocks. After weaning they received soy crop and grain sorgum supplementation until slaughter, with final live weights ranging from 67 to 72 kg at 14 months old. Samples were collected immediately (within 30 minutes) after slaughter from the mid superficial belly of the following muscles: Semimembranosus (SM), Semitendinosus (ST), Cranial Gluteobiceps (GB), Gluteus medius (GM), Rectus femoris (RF), Supraspinatus (SS), Longissimus lumborum (LL), Adductor (AD) and Psoas major (PM). The anatomical location of the muscles sampled is presented as supplementary material (S1 Fig). Samples were embedded in RNAlater solution (Ambion) and stored at -20°C until RNA extraction.
Histochemical tests for fiber typing, fiber morphometric analysis and meat quality determination were performed by Ithurralde et al. (2015, 2017)[12,13].
RNA extraction, library generation and sequencing
Muscle samples were frozen in liquid nitrogen and thoroughly ground to powder using mortar and pestle. Total cellular RNA was isolated and purified with TRIzol® Reagent in conjunction with PureLink RNA Mini Kit (Cat.12183018A). To avoid genomic DNA contamination, eluted RNA was treated with DNase-free DNase Set (QIAGEN). RNA yield and quality were measured through an Agilent 2100 Bioanalyzer (Model G2939B).
Libraries from individual samples were sequenced with Illumina’s HiSeq 2000 at Macrogen (Korea), by means of 100 pb paired-end reads, in 4 lanes and 9 samples per lane. Libraries were barcoded and multiplexed before sequencing.
Mapping reads to the sheep reference genome
Sequencing reads were mapped to the sheep (Ovis aries) reference genome (Oar_v3.1) using the software Tophat (v2.0.13) [14, 15]. Following [16] and [17], two consecutive rounds of alignments were performed in order to maximize sensitivity to splice junction discovery. Briefly, novel splice junctions were initially discovered in each sample, and then these novel splice junctions plus known splice junctions from Ensembl annotation were combined and supplied to Tophat for a second round of alignment. This strategy allows a full utilization of the novel junctions identified in the samples.
Assembly of the transcripts and estimation of the abundance
Transcript models were reconstructed using Cufflinks (v2.2.1) [18]. This software recreates a set of transcript models that best explain the sequencing alignments observed in the samples. The reference ovine annotation (GTF file) was merged with sample assemblies with the purpose of combining annotated (known) transcripts with novel transcript units. This approach maximizes the quality of the final assembly [16, 17]. Finally, in each muscle sample, the total number of reads that effectively mapped to each gene described in the final assembly was calculated using the python script htseq-count [19].
Gene expression analysis
Differentially expressed genes (DEG) between muscles were detected using the R package edgeR (v.3.12.0) [20]. Briefly, edgeR uses the trimmed mean of M-values (TMM) as RNA-Seq normalization method, calculates gene wise negative binomial dispersion parameters using an empirical Bayes procedure, and fits generalized linear models followed by likelihood ratio tests for detecting differentially expressed genes [21]. This study consisted in a paired design because pairwise comparisons were performed between muscles, and each animal in the study posed both muscles. In this context, the likelihood ratio test detected genes that were differentially expressed between muscles, adjusting simultaneously for differences between animals. This test can be viewed as a generalization of a paired t-test. This experimental design increases the statistical power and reduces potential confounder effects. The P-values were adjusted for multiple testing using the Benjamini-Hochberg procedure [22].
Gene and protein functions were taken mainly from Ensembl browser and UniProt database (http://www.uniprot.org/).
Gene set enrichment analysis
The significant enrichment of Gene Ontology (GO) functional terms with genes that showed differential expression between muscles of interest was analyzed using a Fisher’s exact test. Significant genes with FDR ≤ 0.05 and Ensembl annotations were tested against the background set of all expressed genes with Ensembl annotations. The analysis was performed using the R package goseq [23].
Validation by RT-PCR
Eight genes that showed highly significant differential expression between ST and SS (the most contrasting muscles in terms of metabolic and histochemical properties according to Ithurralde et al., 2015, 2017 [12, 13]) were chosen for validation of the RNA-Seq results: myozenin 2 (MYOZ2), myosin light chain 2 (MYL2), myosin light chain 3 (MYL3), tropomyosin 3 (TPM3), myosin light chain 6 (MYL6), troponin T1 (TNNT1), troponin C1 (TNNC1) and methyl transferase-like 21C (METTL21C). The same muscle samples used for RNA-Seq were used for validation of differential expression. The validation was performed using quantitative real-time PCR (qPCR) conducted with theCorbett Rotor-Gene 6000 system (QIAGEN). A total of 1μg RNA from each sample was used to synthesize cDNA using the SuperScript®III Transcriptase (Invitrogen) with primers Oligo-dT (Invitrogen) following the manufacturer’s instructions. Primers (S1 Table) to specifically amplify cDNA of target genes were designed using the Primer3 software (http://bioinfo.ut.ee/primer3/) based on ovine genome sequences (Oar_v3.1/oviAri3) available from UCSC Genome Browser (https://genome.ucsc.edu). Primers were designed to cross exon-exon junctions to minimize the potential of amplifying genomic DNA. Real-time PCR reactions were performed in a total volume of 15 μl containing 2 uL (20 ng) cDNA, 1uL forward and reverse primer mix (300nM each), 4.5 uL PCR water and 7.5 μL QuantiFast SYBR Green 2x qPCR Master Mix, QIAGEN), using the following standard amplification conditions for all sets of primers: 5 min at 95°C and 40 cycles of 15 s at 95°C, 45 s at 60°C and 20 s at 72°C. Dissociation curves were run on all samples to detect primer-dimers, contamination or presence of other amplicons. For each gene assessed, samples were run in triplicate. Each run disc included a triplicate pool of cDNA from all ST and SS samples assayed (exogenous control) and a duplicate non-template control (NTC). Gene ACTB (β-actin) [24] was used as the housekeeping gene (endogenous control). Gene expression was measured by relative quantification (2-ΔΔCt method, [25]) to the exogenous control and normalized to the endogenous control. Amplification efficiencies (E = [10-1/slope-1]) for all genes were estimated by linear regression (Rotor Gene 6000 Q Software 1.7.75) of a triplicate cDNA standard curve (n = 5 dilutions, from 100 to 6.25ng/tube). Amplification efficiencies ranged from 0.94 to 1.10. The intra-assay %CVs between triplicates were < 2%.
Results and discussion
Overall evaluation of gene expression in different muscles
Whole transcriptomes from nine different muscles of four unrelated cross-bred lambs were successfully sequenced using RNA-Seq technology. Sequencing data can be accessed through NCBI GEO with accession number GSE112500. A summary of sequencing read alignments to the reference genome is presented in S2 Table. The number of uniquely mapped reads that were effectively used for the gene expression analysis ranged between 4 and 6 million per sample. This number of mapped reads is enough to accurately quantify moderate to highly expressed genes (e.g., [26]). Indeed, for gene expression, sequencing depth is much less important than replication [27], and so, we preferred to include more animals (i.e. more biological replicates per muscle) than increasing read coverage.
Controlling false discovery rate (FDR) at 1%, Fig 1 shows the total number of genes that showed at least a 2-fold expression difference (Fold Change > 2) in each pairwise comparison. Significant differences were detected between almost all muscle comparisons. While some muscles show very similar expression patterns, with no genes differentially expressed (AD vs. GM) or very few (e.g. 4 genes between GB and GM, 11 genes between AD and LL, and 12 between AD and PM), some others show striking differences; for example, there were 479 differentially expressed genes between SM and SS, 395 between ST and GM, 385 between SS and ST, and a total of 379 between SS and LL.
Controlling false discovery rate (FDR) at 1%, the matrix shows the total number of genes that showed at least a 2-fold expression difference (Fold Change > 2) in each pairwise comparison. SM: Semimembranosus, ST: Semitendinosus, GB: cranial Gluteobiceps, GM: Gluteus medius, RF: Rectus femoris, SS: Supraspinatus, LL: Longissimus lumborum, AD: Adductor, PM: Psoas major. (File: Fig 1.tiff).
Alternative multidimensional scaling (MDS) analyses were performed in order to visualize the overall relationship between the samples under study (Fig 2). The MDS plot displaying the 36 samples clearly shows that the greatest variation is explained by the animal effect, but within each animal, the muscles are usually separated following the same pattern, e.g. SS and ST are consistently very different (Fig 2, left). Note that this is a paired design study, and hence, differentially expressed genes were detected adjusting simultaneously for differences between animals. The MDS plot based on the set of DEG shows that both dimensions clearly separate muscles SS, SM, ST, GM and LL, with the other four muscles in intermediate positions (Fig 2, right).
Using all the samples (left) and using the set of differentially expressed genes (right). SM: Semimembranosus, ST: Semitendinosus, GB: cranial Gluteobiceps, GM: Gluteus medius, RF: Rectus femoris, SS: Supraspinatus, LL: Longissimus lumborum, AD: Adductor, PM: Psoas major. (File: Fig 2.tiff).
There is a great agreement between dimension 1 of the MDS plot and the ATPase (contraction speed) and glycolytic activities of these muscles; in particular, SS has the lowest ATPase activities and also the lowest glycogen synthetase, phosphorylase, and phosphofructokinase activities, while SM and ST have the highest ATPase and also the highest phosphorylase, glycogen synthetase, and phophofructokinase activities [28]. Dimension 2 clearly separates SM, LL and GM from ST and SS. According to [28], SM and Longissimus dorsi are considered fast twitch red muscles, having characteristics typical of fast muscles (high ATPase and high glycolytic activity) but also of slow muscles (high oxidative activity). Ithurralde et al. (2015) found LL (lumbar portion of Longissimus dorsi) within the cluster of fast glycolytic muscles, as well as GM, and clustered SM within fast red muscles [12]. However, the three muscles showed high oxidative fiber proportion and high oxidative scores, as well as high glycolytic fiber proportion. Although being a fast-glycolytic muscle, ST showed extreme values and significant differences in terms of fiber type and morphology in comparison with both GM and LL [12].
Interestingly, MDS plots showed a very similar pattern to the principal component analysis presented by Ithurralde et al. (2017), in which they study the general relationships between meat quality traits and fiber typing across muscles [13]. In agreement with the present gene expression analysis, the muscles clustered in three main categories, with SS, ST and SM clearly separated. SS showed typical values of slow oxidative muscles, SM is an intermediate fast red muscle, and ST is a typical fast glycolytic muscle [12]. As will be seen later on, the observed differences in gene expression profiles could well be related to these metabolic and histochemical properties.
Functional characterization of differently expressed genes between pairs of muscles
The following analyses will focus on the muscles that showed the highest differences in terms of gene expression (ST, SS, SM, LL and GM).
Table 1 shows the 30 most significant differentially expressed genes (threshold: FDR < 1%) between ST and SS, all of them with fold change > 2 (log2FC >1). Several genes that were up regulated in SS are related to muscle tissue structure, including myosin heavy chain 7B and troponin T1, C1 and I1, all common isoforms of slow muscle cells [29], and in agreement with previous transcriptome studies in sheep muscle [6] and in mice [11].
Positive log2 fold change (Log2FC) means higher expression in ST vs. SS.
ATP2A2 and CASQ2 genes are related to calcium ion homeostasis. ATP2A2 enzyme is a slow muscle isoform that transports calcium between the cytosol and the sarcoplasmic reticulum lumen, and is involved in the regulation of contraction and relaxation [30]. CASQ2 is the slow skeletal muscle isoform of calsequestrin, a calcium-binding protein that works as a calcium store in muscle cells, regulating the release of calcium through ryanodine channels [31]. Even if calcium channels are more densely distributed in fast fibers than in slow ones (to accelerate charge movement), calcium entry through the sarcolemma is higher in slow fibers, at rest and during depolarization [29]. Our findings support the importance of calcium homeostasis in a slow twitch muscle as SS, and later on its implications on water holding capacity will be discussed.
Troponins, myosin and ATP2A2 are also involved in the transition between fast and slow fibers [29]. This could be related to muscle plasticity, a property of the muscle fibers that can change their characteristics in response to external influences, specially to nerve activity [29]. Myosin gene expression is also related to fiber type changes and abundance [8]. Up-regulation of genes related to myofiber type and composition in SS shows a higher fiber heterogeneity of slow muscles with respect to fast muscles, as was observed in RNA-Seq studies in mice [11].
On the other hand, genes more expressed in ST are related to morphogenesis (HOXD8 and HOXC10, [32], carbohydrate catabolic processes and protein phosphorylation (CHI3L1, [33], and glucose homeostasis (STXBP5L, [34], among others, all functions associated with high metabolic activity and tissue development in a typical fast-twitch muscle. High expression of genes related to morphogenesis and development were also found in the transcriptome of Chinese lamb muscles [6], and may be related to growth in young animals.
Accordingly, Ithurralde et al. (2015) found, in the histochemical analyses performed in the same muscle samples used for the present study, a typical slow muscle fiber typing profile for SS (the highest percentage of type I fibers and the lowest percentage of type II and glycolytic fibers, combined with a high oxidative score) and a typical fast glycolytic profile for ST (the lowest percentage of type I fibers and of oxidative fibers, but the highest percentage of type II and glycolytic ones) [12].
Two Novel Transcript Units (NTU, unknown/unannotated ovine genes) were detected in ST vs. SS comparison (Table 1). The first NTU of the table (log2FC = -4.38, P value = 3.33e-15) mapped to OAR1 (chr1: 230473614–230478546). This is a conserved region according to Ensembl browser comparative genomics chart, flanked by PLCH1 and MME genes, as well as by two ncRNA. A CDS (Ensembl annotation RP11-451G4.2) and a long intervening no coding RNA (lincRNA; RP11-451G4.3) were reported for this syntenic region in the human genome. According to GO annotations in Ensembl browser and UniProt database, there is RNASeq evidence that the human CDS is differentially expressed in muscle, related to slow-twitch skeletal muscle fiber contraction regulation and calcium ion transport, two functions associated with its up regulation in SS. The other NTU (ESRRG) mapped to OAR 12 (chr12: 18687137–18712068), in a conserved region flanked by ESRRG and GPATCH2 genes. In the human homologous region there is an intron shared by two ESRRG isoforms (Estrogen-Related Receptor γ). This is a remarkable gene as it is considered a key regulator of oxidative metabolism in skeletal muscle, increasing mitochondrial activity. It is up regulated in slow-twitch type I fibers, required for long term adaptation to exercise, and induces expression of several other genes related to calcium ion transport, muscle contraction, fatty acid metabolism and angiogenesis [35]. In the present study we found several genes up regulated in SS directly related to these functions, as can be seen in Table 1.
The significant enrichment of Gene Ontology (GO) terms with differentially expressed genes was evaluated using a Fisher´s exact test. Fig 3 displays the most relevant functional categories that were over-represented with significant genes detected in ST vs. SS pairwise comparison. Important differences were observed in functional terms closely related to energy generation and carbohydrate metabolic processes; as expected, all these terms were up-regulated in ST compared to SS, in accordance with fast-glycolytic metabolism. ST also showed higher expression of genes related to carboxylic acids and lipid metabolism, which may include amino acids and fatty acids related to protein turnover and energy expenditure. As seen from the most DE genes analysis (Table 1), processes related to muscle function and its regulation were also active in both muscles. ST showed up regulation of genes related to cell communication, which includes cell-cell adhesion and attachment to the extracellular matrix, and to biosynthetic processes. Although there were more genes related to calcium transportation within the top 30 DEG in SS than in ST (Table 1), the latter showed more genes associated with response to calcium ions, consistent with the fact that fast glycolytic muscles need to ensure quick charge movements to allow high speed of contraction [29].
The graphs show the number of differentially expressed genes in each functional category. P values are in parentheses. Genes up-regulated in ST are in light red while genes up-regulated in SS are in dark red. (File: Fig 3.tif).
Table 2 shows the most differently expressed genes between SM (fast red muscle) and SS (slow red muscle). Some of the genes that were up regulated in SM are related to morphogenesis and development (HOXD8 and HOXC10 [32]), carbohydrate catabolic processes (CHI3L1, [33]), protein phosphorylation (CHI3L1 and KNDC1 [33, 36]), lactate transport (SLC16A3, [37]), glycogen metabolism (PPP1R3D, [38]) and muscle contraction (SNTB1, [39]); all of them typical features of fast muscles. On the other hand, SS showed higher expression of two genes related to calcium transport (ATP2A2 [30] and PVALB [40]), similar to what was observed in the comparison ST vs. SS, one to sodium transport (ADRB2 [41]), and a gene that encodes a potassium channel, KCNMA1 [42]. Potassium is less important than other ions for muscle contraction, but slow muscle fibers are more active than fast ones in removing potassium from the extracellular space [29]. Both SM and SS showed up regulation of genes related to muscle contraction and cell proliferation. This may be associated with the intermediate characteristics of SM compared to SS and ST, combining contractile and metabolic properties of both fast and slow muscles. In the former fiber typing study by Ithurralde et al. (2015), SM had similar values than ST for type I and type II fibers, but is more similar to SS in terms of oxidative and glycolytic fibers percentage and oxidative score [12].
Positive log2 fold change (Log2FC) means higher expression in SM vs. SS.
There were nine NTUs detected in this pairwise comparison, all except one up regulated in SS (Table 2). Six could be assigned to a particular gene (MYBPC1, slow muscle type myosin binding protein C; ADRB2, beta-2 adrenergic receptor; MYLK4, myosin light chain kinase 4; and RORA, RAR related orphan receptor A). The fourth NTU (the only one that is up-regulated in SM) mapped to an intergenic region near a gene that encodes a collagen chain (COL2A1). The sixth NTU (Log2FC = -1.60; P value = 2.17e-10) mapped to OAR 2 (chr2: 220339844–220346024) in a region flanked by DES gene (desmin, a muscle-specific protein of the sarcomere [43]) and a non-characterized gene (ENSOARG00000020189) that in cattle and human is annotated as SPEG gene, related to protein phosphorylation, muscle cell differentiation and development [44]. The seventh NTU mapped to an uncharacterized region with no annotations nor chromosome assignment in OAR v.3.1. Three NTUs mapped to very close regions on OAR 7 that in the human genome fall within RORA gene. This gene encodes a transcription factor that plays an important role in lipid metabolism, angiogenesis and muscle cell differentiation, especially in red muscles [45].
It is interesting to note that several genes related to the transportation of calcium and other ions were up regulated in SS with respect to SM. This could be related to the significant differences in water holding capacity that were found between these muscles by Ithurralde et al. (2017)[13]. Cooking loss was 34.6% in SS and 27.6% in SM (p value < 0.05) while expressed juice (water expelled after compression) was 15.9% and 10.7%, respectively (p value < 0.05). Accelerated calcium release into the sarcoplasm increases the rate of pH decline, favors protein denaturation and induces low water holding capacity, producing meat of lower quality [46]. Although Ithurralde et al. (2017) found lower final pH values in SM than in SS (p < 0.05), this process is more affected by the rate of decline than by the final value itself [13].
Similar results can be seen for ST vs. SS comparison, were again there were more calcium transport related genes up regulated in SS. One of these genes is CASQ2, a calcium release regulator through ryanodine channels that is associated with PSE (pale, soft and exudative) meat in pigs [31, 46]. Significant differences between SS and ST were found for expressed juice, but not for cooking loss [13]. These genes could be also associated to sarcomere length, as calcium transport to and from the sarcoplasmic reticulum in post-mortem muscle is one of the proposed mechanisms to explain observed differences in muscle susceptibility to cold-shortening [47]. In accordance with this, Ithurralde et al. (2017) found significant differences in sarcomere length between SS and ST. These gene expression patterns could also affect calpains, calcium-dependent proteases responsible for meat tenderization, but no significant differences in Warner-Bratzler shear force (WBSF) were found between SS, SM and ST (p > 0.05; [13]).
Few GO terms were significantly enriched with DEG from the SM vs. SS comparison (Fig 4), in accordance with some similar metabolic and histochemical properties of these muscles [12]. SM showed a very strong up regulation of amide and peptide metabolic processes, especially related to translation and peptide synthesis. These differences could be related to a higher protein turnover and metabolic rate in a fast-red muscle in comparison to a slow-red one. There were also significant genes associated with the circulatory system, probably related to the high vascularization of red muscles [12].
The graphs show the number of differentially expressed genes in each functional category. P values are in parentheses. Genes up-regulated in SM are in light red while genes up-regulated in SS are in dark red (File: Fig 4.tif).
The most significant differentially expressed genes between LL and SS are shown in S3 Table. LL showed some similarities with ST, with respect to the high expression of HOXD8, HOXC10 and CHI3L1, and also to SM with up-regulation of SIM1, SNTB1 and SMYD2, all development-related genes [48, 39, 49]. Being the lumbar portion of Longissimus dorsi muscle, LL can either be considered as a fast-glycolytic muscle [12] or a fast-oxidative muscle [28], which may explain our findings. LL showed higher expression of a potassium channel subunit (KCNS3) that combined with KCNB1, another channel subunit more abundant in slow and fast red fibers, generates an important Na-K pump during repetitive electrical activity [29]. SS showed higher expression of genes related to muscle tissue structure and development, including MYL2, MYL6, MYOZ, PDLIM1 and EYA1 [29, 50]. Two genes related to synaptic transmission were found differentially expressed in SS vs. LL: genes PLCL2 (higher expression in LL) and DLG2 (higher expression in SS) [51, 52]. These genes could be involved in the neuromuscular junction, where motoneurons and muscle fibers interact to make possible muscle contraction and movement, both in slow and fast muscles [29].
Most of the differentially expressed genes in the LL vs. SS comparison were related to metabolic properties (Fig 5). Genes related to carbohydrate metabolism and energy generation were up regulated in LL, as expected in a fast-glycolytic muscle.
The graphs show the number of differentially expressed genes in each functional category. P values are in parentheses. Genes up-regulated in LL are in light red while genes up-regulated in SS are in dark red. (File: Fig 5.tif).
With regard to the comparison ST vs. GM (S4 Table), several up regulated genes in ST are consistent with fast white muscles characteristics, including muscle contraction (LMOD1, [53] and MYLK4 [54]), development (PITX2 [55]), glucose homeostasis (STXBP5L [34]), signal transduction (MPP6 [56]), protein phosphorylation and energy metabolism (PRKAB2 [57]). Some of the genes up regulated in GM resemble the expression pattern of red muscles, such as MYH7B, MYOZ2, RRAD (the latter related to glycolysis inhibition and calcium channels regulation [58]) and VEGFA (vascular growth factor [59]. Most GO expression differences in ST vs. GM comparison were not associated with energy metabolism but to catalytic activity regulation and response to external stimuli (Fig 6).
The graphs show the number of differentially expressed genes in each functional category. P values are in parentheses. Genes up-regulated in ST are in light red while genes up-regulated in GM are in dark red. (File: Fig 6.tif).
With regard to meat quality, Ithurralde et al. (2017) reported that redness parameter a* was significantly higher (p < 0.05) in SS than in ST, LL, and GM, probably related to its oxidative metabolism [13]. SS also showed the highest yellowness parameter b*, followed by ST, which the authors relate to higher amounts of intramuscular fat and connective tissue with respect to other muscles. In this respect, it is interesting that FABP3, a gene that encodes a fatty acid transport enzyme, is one of the top 30 differentially expressed genes in SS (Table 1) and that lipid metabolic processes were up-regulated in ST (Fig 3).
Several genes that were more expressed in SS are related to cell adhesion and cytoskeleton organization (PDLIM1, CCDC63, TPM3, MAST2 and TACC2 [60, 29, 61, 62]) and extracellular matrix components (FREM2 [63]), which in turn are associated with the high connective tissue content, high yellowness parameter and relatively high WBSF score (2.28 kg/cm2) found by Ithurralde et al. (2017) for this muscle [13]). Accordingly, in a comparison between 18 lamb muscles, [64] found that SS was one of the muscles with higher collagen content. The relatively high expression of genes involved in muscle contraction and development in SS (Figs 3–5) could be related to this.
Protein phosphorylation is negatively related to post-mortem tenderization [65]. Several genes involved in protein phosphorylation were found up regulated in SS, SM and ST, three muscles that showed high WBSF scores (2.28, 2.73 and 2.55 kg/cm2, respectively), in contrast with GM, a muscle with low WBSF (1.81 kg/cm2; [12]).
Higher glycolytic activities and lactic acid generation of fast or intermediate muscles tend to decrease meat pH values, being significantly lower in ST, LL, SM and GM than in SS (p < 0.05; [12]). This is in agreement with our findings that several carbohydrate and carboxylic acids (e.g. lactic acid) metabolic processes were up regulated in ST in comparison with SS (Fig 3), and with the high expression of SLC16A3 in SM, a gene closely involved in lactate transmembrane transport [37].
Validation of gene expression analysis
In order to validate the findings of the RNASeq analysis, the expression of eight genes, namely MYOZ2, MYL2, MYL3, TPM3, MYL6, TNNT1, TNNC1 and METTL21C, was assessed using qRT-PCR. Fig 7 shows the fold differences in gene expression measured by both RNA-Seq and qRT-PCR, confirming that the eight genes showed similar patterns of mRNA abundance with both methods.
The eight genes showed higher expression in Supraspinatus (SS) vs. Semitendinosus (ST). Spearman correlation value = 0.84.
Conclusions
Gene expression profiles obtained with RNA-Seq technology clearly separated muscles with different histochemical and metabolic properties. Pairwise comparisons between muscles show significant differences in terms of both major genes and relevant functions and processes closely related to meat quality parameters. Our results reveal new candidate genes associated with meat quality in sheep that deserve future research, and give a deeper insight into the genetic architecture of these complex traits.
Supporting information
S1 Fig. Schematic diagram showing the anatomical position of the studied muscles in a lamb carcass: SM: muscle Semimembranosus; GB: muscle Gluteobiceps; ST: muscle Semitendinosus; GM, muscle Gluteus medius; SS: muscle Supraspinatus.
a) Transverse section of the hind limb showing muscles Rectus femoris (RF, cranial position), Gluteobiceps(GB, lateral position), Semitendinosus (ST, caudo-lateral position), Semimembranosus (SM caudo-medial position) and Adductor (AD, medial position). b) Transverse section of the loin showing muscles Longissimus lumborum (LL, dorsal position) and Psoas major (PM, ventral position).
https://doi.org/10.1371/journal.pone.0200732.s001
(DOCX)
S1 Table. Primers used for qPCR validation analysis.
https://doi.org/10.1371/journal.pone.0200732.s002
(DOCX)
S2 Table. Summary of sequence read alignments to the reference genome.
https://doi.org/10.1371/journal.pone.0200732.s003
(XLSX)
S3 Table. Top 30 most significant genes detected in Longissimus lumborum (LL) vs. Supraspinatus (SS) pairwise comparison.
Positive log2 fold change (Log2FC) means higher expression in LM vs. SS.
https://doi.org/10.1371/journal.pone.0200732.s004
(DOCX)
S4 Table. Genes contributing to the Gene Ontology terms analyzed.
Gene name: gene Ensembl ID.
https://doi.org/10.1371/journal.pone.0200732.s005
(XLSX)
Acknowledgments
Computational resources were provided by the Florida Agricultural Experiment Station and the Department of Animal Sciences, University of Florida.
References
- 1. Knight MI, Daetwyler HD, Hayes BJ, Hayden MJ, Ball AJ, Pethick DW, et al. An independent validation association study of carcass quality, shear force, intramuscular fat percentage and omega-3 polyunsaturated fatty acid content with gene markers in Australian lamb. Meat Sci. 2014 Feb; 96(2 Pt B):1025–33. Epub 2013 Jul 17. pmid:23948658
- 2. Bernard C, Cassar-Malek I, Le Cunff M, Dubroeucq H, Renand G, Hocquette JF. New indicators of beef sensory quality revealed by expression of specific genes. J Agric Food Chem. 2007 Jun 27;55(13):5229–37. Epub 2007 Jun 5. pmid:17547415
- 3. Hocquette JF, Botreau R, Picard B, Jacquet A, Pethick DW, Scollan ND. Opportunities for predicting and manipulating beef quality. Meat Sci. 2012 Nov;92(3):197–209. Epub 2012 Apr 12. pmid:22554471
- 4. Damon M, Denieul K, Vincent A, Bonhomme N, Wyszynska-Koko J, Lebret B. Associations between muscle gene expression pattern and technological and sensory meat traits highlight new biomarkers for pork quality assessment. Meat Sci. 2013 Nov;95(3):744–54. Epub 2013 Feb 15. pmid:23481319
- 5. Zhang C, Wang G, Wang J, Ji Z, Liu Z, Pi X, et al. Characterization and comparative analyses of muscle transcriptomes in Dorper and Small-Tailed Han sheep using RNA-Seq technique. PLoS One. 2013 Aug 30;8(8):e72686. eCollection 2013. pmid:24023632
- 6. Zhang C, Wang G, Wang J, Ji Z, Dong F, Chao T. Analysis of differential gene expression and novel transcript units of ovine muscle transcriptomes. PLoS One. 2014 Feb 26;9(2):e89817. eCollection 2014. pmid:24587058
- 7. Sun L, Bai M, Xiang L, Zhang G, Ma W, Jiang H. Comparative transcriptome profiling of longissimus muscle tissues from Qianhua Mutton Merino and Small Tail Han sheep. Sci Rep. 2016 Sep 20;6:33586. pmid:27645777
- 8. Vuocolo T, Byrne K, White J, McWilliam S, Reverter A, Cockett NE, et al. Identification of a gene network contributing to hypertrophy in callipyge skeletal muscle. Physiol Genomics. 2007 Feb 12;28(3):253–72. Epub 2006 Oct 31. pmid:17077277
- 9. Fleming-Waddell JN, Olbricht GR, Taxis TM, White JD, Vuocolo T, Craig B A, et al. Effect of DLK1 and RTL1 but not MEG3 or MEG8 on muscle gene expression in Callipyge lambs. PLoS One. 2009 Oct 9;4(10):e7399. pmid:19816583
- 10. Zhang W, Wang L, Zhou P, Song G, Shen M, Gan S, et al. Identification and analysis of genetic variations in pri-miRNAs expressed specifically or at a high level in sheep skeletal muscle. PLoS One. 2015 Feb 20;10(2):e0117327. eCollection 2015. pmid:25699993
- 11. Raz V, Riaz M, Tatum Z, Kielbasa SM, 't Hoen PAC. The distinct transcriptomes of slow and fast adult muscles are delineated by noncoding RNAs. FASEB J. 2017 Nov 15. pii: fj.201700861R. pmid:29141996
- 12. Ithurralde J, Bianchi G, Feed O, Nan F, Garibotto G, Bielli A. Histochemical fiber types in 16 heavy-lamb skeletal muscles. Small Rumin Res. 2015 Apr 125: 88–92.
- 13. Ithurralde J, Bianchi G, Feed O, Nan F, Ballesteros F, Garibotto G, et al. Variation in instrumental meat quality among 15 muscles from 14-month-old sheep and its relationship with fibre typing. Anim. Prod. Sci. 2017.
- 14. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009 May 1;25(9):1105–11. Epub 2009 Mar 16. pmid:19289445
- 15. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013 Apr 25;14(4):R36. pmid:23618408
- 16. Peñagaricano F, Souza AH, Carvalho PD, Driver AM, Gambra R, Kropp J, et al. Effect of maternal methionine supplementation on the transcriptome of bovine preimplantation embryos. PLoS One. 2013 Aug 21;8(8):e72302. eCollection 2013. pmid:23991086
- 17. Peñagaricano F, Wang X, Rosa G, Radunz A, Khatib H. Maternal nutrition induces gene expression changes in fetal muscle and adipose tissues in sheep. BMC Genomics. 2014 Nov 28;15:1034. pmid:25429728
- 18. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010 May;28(5):511–5. Epub 2010 May 2. pmid:20436464
- 19. Anders SP, Pyl T, Huber W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015 Jan 15;31(2):166–9. Epub 2014 Sep 25. pmid:25260700
- 20. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010 Jan 1;26(1):139–40. Epub 2009 Nov 11. pmid:19910308
- 21. McCarthy DJ, Chen YS, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012 May;40(10):4288–97. Epub 2012 Jan 28. pmid:22287627
- 22. Benjamini Y, Hochberg Y. Controlling the false discovery rate—a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B. 1995 57: 289–300.
- 23. Young M, Wakefield M, Smyth G, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14. Epub 2010 Feb 4. pmid:20132535
- 24. Vorachek WR, Hugejiletu H, Bobe G, Hall JA. Reference gene selection for quantitative PCR studies in sheep neutrophils. Int J Mol Sci. 2013 May 30;14(6):11484–95. pmid:23722658
- 25. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001 May 1;29(9):e45. pmid:11328886
- 26. Lei R, Ye K, Gu Z, Sun X. Diminishing returns in next-generation sequencing (NGS) transcriptome data. Gene. 2015 Feb 15;557(1):82–7. Epub 2014 Dec 10. pmid:25497830
- 27. Liu Y, Zhou J and White K.P. RNA-seq differential expression studies: more sequence or more replication? Bioinformatics. 2014 Feb 1; 30(3): 301–304. pmid:24319002
- 28. Briand M, Talmant A, Briand Y, Monin G, Durand R. Metabolic Types of Muscle in the Sheep: I. Myosin ATPase, Glycolytic, and Mitochondrial Enzyme Activities. Eur J Appl Physiol Occup Physiol. 1981;46(4):347–58. pmid:6455290
- 29. Schiaffino S, Reggiani C. Fiber types in mammalian skeletal muscles. Physiol Rev. 2011 Oct;91(4):1447–531. pmid:22013216
- 30. Dally S, Bredoux R, Corvazier E, Andersen JP, Clausen JD, Dode L, et al. Ca2+-ATPases in non-failing and failing heart: evidence for a novel cardiac sarco/endoplasmic reticulum Ca2+-ATPase 2 isoform (SERCA2c). Biochem. J. 2006; 395:249–258. pmid:16402920
- 31. Sanchez EJ, Munske GR, Criswell A, Milting H, Dunker AK, Kang C. Phosphorylation of human calsequestrin: implications for calcium regulation. Mol. Cell. Biochem. 2011; 353:195–204. pmid:21416293
- 32. Murtha MT, Leckman JF, Ruddle FH. Detection of homeobox genes in development and evolution. Proc. Natl. Acad. Sci. U.S.A. 1991; 88:10711–10715. pmid:1720547
- 33. Chen CC, Llado V, Eurich K, Tran HT, Mizoguchi E. Carbohydrate-binding motif in chitinase 3-like 1 (CHI3L1/YKL-40) specifically activates Akt signaling pathway in colonic epithelial cells. Clin. Immunol. 2011; 140:268–275. pmid:21546314
- 34. Bhatnagar S, Oler AT, Rabaglia ME, Stapleton DS, Schueler KL, Truchan NA, et al. Positional cloning of a type 2 diabetes quantitative trait locus; tomosyn-2, a negative regulator of insulin secretion. PLoS Genet. 2011; 7:e1002323–e1002323. pmid:21998599
- 35. Rangwala SM, Wang X, Calvo JA, Lindsley L, Zhang Y, Deyneko G, et al. Estrogen-related receptor is a key regulator of muscle mitochondrial activity and oxidative capacity. J Biol Chem. 2010 Jul 16;285(29):22619–29. Epub 2010 Apr 23. pmid:20418374
- 36. Huang J, Furuya A, Furuichi T. Very-KIND, a KIND domain containing RasGEF, controls dendrite growth by linking Ras small GTPases and MAP2. J. Cell Biol. 2007; 179:539–552. pmid:17984326
- 37. Wilson MC, Meredith D, Fox JE, Manoharan C, Davies AJ, Halestrap AP. Basigin (CD147) is the target for organomercurial inhibition of monocarboxylate transporter isoforms 1 and 4: the ancillary protein for the insensitive MCT2 is EMBIGIN (gp70). J. Biol. Chem. 2005; 280:27213–27221. pmid:15917240
- 38. Armstrong C.G., Browne G.J., Cohen P., Cohen P.T.W. PPP1R6, a novel member of the family of glycogen-targeting subunits of protein phosphatase 1. FEBS Lett. 1997; 418:210–214. pmid:9414128
- 39. Ahn AH, Yoshida M, Anderson MS, Feener CA, Selig S, Hagiwara Y, et al. Cloning of human basic A1, a distinct 59-kDa dystrophin-associated protein encoded on chromosome 8q23-24. Proc. Natl. Acad. Sci. U.S.A. 1994; 91:4446–4450. pmid:8183929
- 40. Troxler H, Kuster T, Rhyner JA, Gehrig P, Heizmann CW. Electrospray ionization mass spectrometry: analysis of the Ca2+-binding properties of human recombinant alpha-parvalbumin and nine mutant proteins. Anal. Biochem. 1999; 268:64–71. pmid:10036163
- 41. Mutlu GM, Dumasius V, Burhop J, McShane PJ, Meng FJ, Welch L, et al. Upregulation of alveolar epithelial active Na+ transport is dependent on beta2-adrenergic receptor signaling. Circ. Res. 2004; 94:1091–1100. pmid:15016730
- 42. Krick S, Platoshyn O, Sweeney M, Kim H, Yuan JX. Activation of K+ channels induces apoptosis in vascular smooth muscle cells. Am. J. Physiol. Cell Physiol. 2001; Apr;280(4):C970–9. pmid:11245614
- 43. Hnia K, Ramspacher C, Vermot J, Laporte J. Desmin in muscle and associated diseases: beyond the structural function. Cell Tissue Res. 2015; 360:591–608. pmid:25358400
- 44. Hsieh CM, Yoshizumi M, Endege WO, Kho CJ, Jain MK, Kashiki S, et al. APEG-1, a novel gene preferentially expressed in aortic smooth muscle cells, is down-regulated by vascular injury. Biol. Chem. 1996; 271:17354–17359.
- 45. Sevane N, Armstrong E, Wiener P, Pong Wong R, Dunner S, GemQual Consortium. Polymorphisms in twelve candidate genes are associated with growth, muscle lipid profile and meat quality traits in eleven European cattle breeds. Mol Biol Rep. 2014 Jul;41(7):4721–31. Epub 2014 Apr 10. pmid:24718780
- 46. Huff-Lonergan E, Lonergan SM. Mechanisms of water-holding capacity of meat: The role of post-mortem biochemical and structural changes. Meat Sci. 2005 Sep;71(1):194–204. Epub 2005 Jun 14. pmid:22064064
- 47. Cornforth DP, Pearson AM, Merkel RA. Relationship of mitochondria and sarcoplasmic reticulum to cold-shortening. Meat Sci. 1980 Apr;4(2):103–21. pmid:22055618
- 48. Coumailleau P, Duprez D. Sim1 and Sim2 expression during chick and mouse limb development. Int. J. Dev. Biol. 2009; 53, 149–157. pmid:19123137
- 49. Diehl F, Brown MA, van Amerongen MJ, Novoyatleva T, Wietelmann A, Harriss J, et al. Cardiac deletion of Smyd2 is dispensable for mouse heart development. PLoS ONE 2010; 5:E9748–E9748. pmid:20305823
- 50. Li X, Oghi KA, Zhang J, Krones A, Bush KT, Glass CK, et al. Eya protein phosphatase activity regulates Six1-Dach-Eya transcriptional effects in mammalian organogenesis. Nature 2003 Nov 20; 426(6964):247–54. pmid:14628042
- 51. Jovanovic JN, Takenaka K, Nakayama KI, Fukami K, Takenawa T, Moss SJ, et al. Modulation of GABA(A) receptor phosphorylation and membrane trafficking by phospholipase C-related inactive protein/protein phosphatase 1 and 2A signaling complex underlying brain-derived neurotrophic factor-dependent regulation of GABAergic inhibition. J. Biol. Chem. 2006; 281:22180–22189. pmid:16754670
- 52. Tao YX, Rumbaugh G, Wang GD, Petralia RS, Zhao C, Kauer FW, et al. Impaired NMDA receptor-mediated postsynaptic function and blunted NMDA receptor-dependent persistent pain in mice lacking postsynaptic density-93 protein. J. Neurosci. 2003; 23:6703–6712. pmid:12890763
- 53. Nanda V, Miano JM. Leiomodin 1, a new serum response factor-dependent target gene expressed preferentially in differentiated smooth muscle cells. J. Biol. Chem. 2012; 287:2459–2467. pmid:22157009
- 54. Birukov KG, Csortos C, Marzilli L, Dudek S, Ma SF, Bresnick AR, et al. Differential regulation of alternatively spliced endothelial cell myosin light chain kinase isoforms by p60(Src). J. Biol. Chem. 2001; 276:8567–8573. pmid:11113114
- 55. Martinez-Fernandez S, Hernandez-Torres F, Franco D, Lyons GE, Navarro F, Aranega AE. Pitx2c overexpression promotes cell proliferation and arrests differentiation in myoblasts. Dev. Dyn. 2006; 235:2930–2939. pmid:16958127
- 56. Leonoudakis D, Conti LR, Anderson S, Radeke CM, McGuire LM, Adams ME, et al. Protein trafficking and anchoring complexes revealed by proteomic analysis of inward rectifier potassium channel (Kir2.x)-associated proteins. J. Biol. Chem. 2004; 279:22331–22346. pmid:15024025
- 57. Thornton C, Snowden MA, Carling D. Identification of a novel AMP-activated protein kinase beta subunit isoform that is highly expressed in skeletal muscle. J. Biol. Chem. 1998; 273:12443–12450. pmid:9575201
- 58. Chang L, Zhang J, Tseng YH, Xie CQ, Ilany J, Bruning JC, et al. Rad GTPase deficiency leads to cardiac hypertrophy. Circulation 2007; 116:2976–2983. pmid:18056528
- 59. Krilleke D, DeErkenez A, Schubert W, Giri I, Robinson GS, Ng YS, et al. Molecular mapping and functional characterization of the VEGF164 heparin-binding domain. J. Biol. Chem. 2007; 282:28045–28056. pmid:17626017
- 60. Kotaka M, Kostin S, Ngai S, Chan K, Lau Y, Lee SM, et al. Interaction of hCLIM1, an enigma family protein, with alpha-actinin 2. J. Cell. Biochem. 2000; 78:558–565. pmid:10861853
- 61. Lumeng C, Phelps S, Crawford GE, Walden PD, Barald K, Chamberlain JS. Interactions between beta 2-syntrophin and a family of microtubule-associated serine/threonine kinases. Nat. Neurosci. 1999; 2:611–617. pmid:10404183
- 62. Dou Z, Ding X, Zereshki A, Zhang Y, Zhang J, Wang F, et al. TTK kinase is essential for the centrosomal localization of TACC2. FEBS Lett. 2004; 572:51–56. pmid:15304323
- 63. Jadeja S, Smyth I, Pitera JE, Taylor MS, van Haelst M, Bentley E, et al. Identification of a new gene mutated in Fraser syndrome and mouse myelencephalic blebs. Nat. Genet. 2005; 37:520–525. pmid:15838507
- 64. Tschirhart-Hoelscher TE, Baird BE, King DA, McKenna DR, Savell JW. Physical, chemical, and histological characteristics of 18 lamb muscles. Meat Sci. 2006 May;73(1):48–54. Epub 2005 Dec 20. pmid:22062053
- 65. Chen L, Xin L, Ni N, Liu Y, Chen L, Wang Z, et al. Phosphorylation of myofibrillar proteins in post-mortem ovine muscle with different tenderness. J Sci Food Agric. 2016 Mar 30;96(5):1474–83. Epub 2015 Jun 15. pmid:25950868