Search tips
Search criteria 


Logo of springeropenLink to Publisher's site
Functional & Integrative Genomics
Funct Integr Genomics. 2010 November; 10(4): 609–618.
Published online 2010 June 4. doi:  10.1007/s10142-010-0175-2
PMCID: PMC2990504

Skeletal muscle specific genes networks in cattle


While physiological differences across skeletal muscles have been described, the differential gene expression underlying them and the discovery of how they interact to perform specific biological processes are largely to be elucidated. The purpose of the present study was, firstly, to profile by cDNA microarrays the differential gene expression between two skeletal muscle types, Psoas major (PM) and Flexor digitorum (FD), in beef cattle and then to interpret the results in the context of a bovine gene coexpression network, detecting possible changes in connectivity across the skeletal muscle system. Eighty four genes were differentially expressed (DE) between muscles. Approximately 54% encoded metabolic enzymes and structural-contractile proteins. DE genes were involved in similar processes and functions, but the proportion of genes in each category varied within each muscle. A correlation matrix was obtained for 61 out of the 84 DE genes from a gene coexpression network. Different groups of coexpression were observed, the largest one having 28 metabolic and contractile genes, up-regulated in PM, and mainly encoding fast-glycolytic fibre structural components and glycolytic enzymes. In FD, genes related to cell support seemed to constitute its identity feature and did not positively correlate to the rest of DE genes in FD. Moreover, changes in connectivity for some DE genes were observed in the different gene ontologies. Our results confirm the existence of a muscle dependent transcription and coexpression pattern and suggest the necessity of integrating different muscle types to perform comprehensive networks for the transcriptional landscape of bovine skeletal muscle.

Electronic supplementary material

The online version of this article (doi:10.1007/s10142-010-0175-2) contains supplementary material, which is available to authorized users.

Keywords: Cattle, Myofibre type, Gene expression, Microarray, Gene coexpression network, Skeletal muscle


The knowledge of how genes regulate differences in skeletal muscle composition and metabolism is a subject of outstanding interest in medicine and agriculture, as the type of muscle plays a role in the development of chronic metabolic diseases (Tanner et al. 2002), as well as in the organoleptic properties of the meat (Thompson 2002). Muscles consist of a heterogeneous population of fibres with different molecular, structural, contractile, and metabolic properties, which contribute to a wide variety of functional capabilities. In the bovine, postnatal muscles have three major fibre types encoded by the expression of three myosin heavy chain (MyHC) genes: the slow-oxidative one or type I, and two fast ones or types IIA and IIX (Maccatrozzo et al. 2004; Tanabe et al. 1998). Slow muscles are mainly composed of type I fibres, oriented to continuous but modest activities mediated by an aerobic metabolism. On the other hand, fast-glycolytic muscles consist of a large proportion of type II fibres and are sporadically demanded for short periods of intense muscular activity associated with anaerobic metabolism (Brent and Tabin 2004). The composition of muscle tissue in bovine is known to be influenced by both the genetic background and the development prior to birth of the individuals (Lehnert et al. 2007). Muscle specialisation is the result of the coordinated expression of contractile and metabolic proteins together with the histological features that characterise the fibres. This coordination has been described in bovines (Moreno-Sánchez et al. 2008) and in some other species (Acevedo and Rivero 2006; Quiroz-Rothe and Rivero 2004). Besides the changes in the expression of the tissue-specific myosin isoforms, skeletal myofibres have the ability to undergo adaptive metabolic changes in response to external stimulus such as nutrition (Cassar-Malek et al. 2004) and exercise (Flück 2006). The identification of differentially expressed (DE) genes involved in the functionality of muscles (contractile, metabolic, and structural properties) and their potential relationship to differences in muscle plasticity and meat quality are recent areas of interest (Bernard et al. 2007).

Microarray technology is a useful tool for the comparative study of gene expression differences and to ascertain those genes involved in the phenotypic characteristics of muscles. In cattle, although the global gene expression differences between muscles has not yet been studied using microarrays, a comparison between Rectus abdominis and Semitendinosus was performed based on human macroarrays (Sudre et al. 2003, 2005). In addition, different microarray approaches have been developed to study the differential gene expression between phenotypically distinct muscles in swine (Bai et al. 2003) and in mice (Campbell et al. 2001; Wu et al. 2003). Genes encoding for structural proteins and metabolic enzymes were the main DE groups in all these experiments. Other groups of DE genes included mitochondrial genes, transcription factors, and genes involved in Ca2+ channels. Moreover, ligands for receptor signals and extracellular matrix related genes were DE in some of these studies, suggesting the presence of species- or muscle-specific drivers of molecular characteristics of muscle.

Once the genes of interest have been detected, the next major challenge of genomics is to elucidate the functions of these genes and to discover how they interact to perform specific biological processes (Stuart et al. 2003). Gene networks are a promising tool for modelling, analysis, and visualisation, being considered a semi-quantitative graphical representation of transcriptional regulation that reveals groups of functionally related genes and their regulators, and genes with high transcriptional connectivity (Hudson et al. 2009). One method for building biological networks is to establish connections between genes whose expression profiles are significantly correlated. Some authors have reported that genes participating in the same pathway are often correlated under a large number of diverse conditions (Reverter et al. 2006a; Segal et al. 2003). In fact, different transcriptional regulatory networks have been proposed attending to both the species-specific coexpression patterns (Gardner et al. 2003) and the across-species ones (Stuart et al. 2003), with several features being species specific. In cattle, two gene coexpression networks for skeletal muscle (from Longissimus dorsi data) have been inferred from microarray gene expression data coming from different experimental conditions using a cDNA platform (Reverter et al. 2006b) or an oligonucleotide one (Hudson et al. 2009). To date, these networks represent the most comprehensive and largest bovine muscle gene profiling data set.

In the literature, the majority of studies on beef cattle concerns Longissimus muscle, but how to extrapolate these results to other muscle types remains challenging. The purpose of the present study was to profile the differential gene expression between two skeletal muscle types in beef cattle, Psoas major (PM) and Flexor digitorum (FD), after slaughter. To this end, three approaches are followed, including: the interpretation of the results in the light of the above-mentioned coexpression network (cDNA platform) of the Longissimus dorsi muscle of cattle, the assessment of the functional relationships among the DE genes and the detection of possible changes in connectivity across the skeletal muscle system. Will the relationships that exist between genes within a single muscle and species be conserved across muscle types?

Materials and methods

Animal material

Ten unrelated male calves of the Avileña-Negra Ibérica beef cattle breed were used for the experiment. The animals were fattened with the same diet during a finishing period of around 200 days and slaughtered when they reached commercial requirements, at an average live weight of 500 kg and 15 months of age. Tissue samples of PM and FD were collected at slaughter, immediately frozen in liquid nitrogen, and stored at −80°C until RNA isolation. Animal Care and Use Committee approval was not obtained for this study because the samples were collected from carcasses. In a previous study, we characterised the fibre composition of PM and FD, in Avileña-Negra Ibérica male calves (Moreno-Sánchez et al. 2008). These two muscles were chosen because of their physiological and meat quality differences. FD exhibited a large proportion of type I (54%) and IIA (36%) fibres and lacked type IIX, while PM was composed by a large proportion of type II (21% IIA and 35% IIX) and a noticeable percentage of type I (29%) fibres.

For validation purposes using real-time RT-PCR, PM and FD samples were obtained from another 15 unrelated male calves, slaughtered under similar conditions as the individuals used in the microarray experiment. Tissue samples of PM and FD were collected and stored as described above.

Microarray methods and data acquisition

The bovine fat and muscle cDNA microarray (Lehnert et al. 2004; Lehnert et al. 2006), arguably the most widely used in cattle muscle gene expression studies, was used in this experiment. While previous studies using the same platform revealed it to be well suited for the expression profiling of skeletal muscle (Byrne et al. 2005; Lehnert et al. 2006, 2007; Reverter et al. 2005, 2006a, 2008; Tan et al. 2006; Wang et al. 2005, 2009), we further addressed the coverage of the transcriptome that could be achieved by the following approach: We downloaded the most recent version (Release 8, July 6, 2009) of the “Flat Files” from the Human Protein Reference Database (HPRD; containing features of proteins such as post-translational modifications, subcellular localisation, protein–protein interactions and tissue expression. Of the 17,618 unique genes represented in the HPRD dataset, 3,925 (22.3%) are reported to be expressed in skeletal muscle. Also, 731 of the 841 genes in our platform are included in the 17,618 HPRD genes. Of these, 231 (31.6%) are reported to be expressed in skeletal muscle. These figures result in an over-representation hypergeometric test p value of 6.97E−10. Although oligonucleotide microarrays are nowadays becoming the standard in human, mouse and in a growing number of species, we concluded that the muscle specificity of the microarray platform used in this study makes it suitable for the expression profiling of bovine skeletal muscle.

Twenty microarray slides were hybridised following a loop design (Kerr 2003) shown in Supplementary Fig. S1. In each slide, cDNA coming from the two types of muscles were simultaneously hybridised muscle samples belonging alternatively to the same or to two different individuals. Dye channel, red (Cy3) or green (Cy5) was swapped between muscles from one slide to the next to account for a possible dye bias.

Total RNA was isolated from frozen muscle tissue using the Qiagen RNeasy Kit (Qiagen, Hilden, Germany). RNA purity and integrity were assessed by RNA Nano Labchips in an Agilent 2100B Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). The cDNA was synthesised from 30 μg of total RNA with the CyScribe™ First-Strand cDNA Labelling Kit (Amersham Biosciences), incorporating either Cy3-dUTP or Cy5-dUTP into the cDNA. Slides were pre-hybridised in a solution containing 6× SSC, 0.5% SDS and 1% bovine serum albumin (Sigma) for 1 h, and then hybridised overnight at 42°C in a solution containing 50% formamide, 6× SSC, 0.5% SDS, 5× Denhardt’s, 20 μg of poly(A) (Sigma) and salmon sperm (Invitrogen). Afterwards, the slides were washed consecutively in three washing solutions (first 0.2× SSC and 0.1% SDS, second 0.2× SSC and third 0.06× SSC) and were finally dried. All the hybridisation steps were carried out in a Lucidea™ SlidePro Automated Hybridization Station (Amersham Biosciences).

Microarray slides were scanned with the GenePix 4000B scanner (Axon Instruments, Union City, CA, USA) at a resolution of 5 μm (PMT values ranging from 550 to 700 and laser power 100%). The intensity value for the Cy3 and Cy5 channels from each spot was acquired by the Gene Pix Pro 5.0 software (Axon Instruments). Invalid spots were filtered out and not included in the further normalisation and analysis processes (Díaz et al. 2009). A total of 269,712 intensity readings from both channels and corresponding to 8,538 unique clones were retained for the analysis. Separate for each channel, intensity readings (y) were background-corrected by subtracting the background (Bg) from the foreground (Fg) intensities and log2-transformed to approximate a normal distribution as follows:

equation M1

Data normalisation and analysis, and clone annotation

Global normalisation and analysis were jointly performed using a Bayesian approach, detailed in a previous study (Díaz et al. 2009). Intensity records were analysed using the model that proved to have the best goodness of fit and the best predicting ability.

A standardised measure of the differences in gene expression for each clone (Δg) was obtained as:

equation M2

where equation M3 and equation M4 represent the posterior means of the effect of the gth clone in the PM and FD muscles, respectively, and sdg represents the standard deviation of the posterior differences.

To determine which clones were DE, a model-based cluster analysis was applied to the Δg values. This analysis grouped clones according to their Δg. The computer package BAYESMIX (Reverter et al. 2003) was used. The program provided posterior probabilities of belonging to each of the clusters for each clone and statistical criteria (AIC, BIC) that allowed the determination of the number of clusters that yielded the best fit to the data. Estimates of false discovery rate (FDR) were obtained from the mixture models as in McLachlan et al. (2005), transforming the n finite numbers of clusters into a two-components mixture, one containing clones that are identified as potentially DE and another group formed by the potentially non-DE clones. Thus,

equation M5

where equation M6 was the posterior probability that the jth gene from the subgroup of DE clones (Nr) was not DE. This has been termed as local false discovery rate (local FDR) (Efron and Tibshirani 2002) and might be seen as an empirical Bayes version of the Benjamini and Hochberg methodology (Efron 2004). The set of DE clones was finally obtained by establishing a desired FDR for the experiment.

Each DE clone was annotated after BLAST searches (Altschul et al. 1990) against different available databases. Gene ontology annotation and enrichment analysis of functional annotation were performed using the David database (Huang et al. 2009) and GOrilla (Eden et al. 2009).

Validation by TaqMan real-time RT-PCR

Real-time RT-PCR, with TaqMan MGB probes, was performed on 11 genes, AK1, ATF4, CA3, CREG1, CRYAB, CSRP3, GPD1, LDHB, MYH1, PDLIM7 and PFKM (see gene names in Table 1). Genes were chosen according to two different criteria, such as to be located in the region corresponding to FDR values between 5% and 10%, or to have a biological relevance for further studies. ACTA1 (Actin alpha1) was used as the internal reference gene as it showed negligible expression differences along the microarray experiment. The sequences of primers and probes, which were designed with the Primer Express software (Applied Biosystems, Foster City, CA, USA), are listed in the Supplementary Table S1. In silico specificity of the amplicons was screened by BLAST searches.

Table 1
Genes up-regulated in PM and FD according to their GO categories. Total number of genes up-regulated in each muscle is shown in brackets. Percentages on the left and right of each GO group are related to the number of genes up-regulated in PM and FD, ...

Real Total ARN Spin Plus kit (Durviz, Valencia, Spain) was used for DNA-free total RNA isolation. RNA quantification and quality control were performed using NanoDrop ND-1000 UV/Vis (Nanodrop Technologies). DNA contamination absence was determined by control PCR in which the reverse transcription step was suppressed. The reactions were set by TaqMan One-Step RT-PCR Master Mix kit (Applied Biosystems). RT-PCR amplification reactions (25 µl) contained 40 ng total RNA, 1× Master Mix without UNG, 1× MultiScribe™ and RNase Inhibitor Mix, 100 nM each of the forward and reverse primers for the internal reference gene, 100 nM each of the forward and reverse primers for the tested gene, and 100 nM each of the TaqMan probes. Reactions were run on an ABI PRISM 7500 Fast Sequence Detector (Applied Biosystems) following the manufacturer’s cycling parameters. Each reaction was repeated three times. The corresponding mRNA levels were measured and analysed by the 7500 System Software (Applied Biosystems). Statistical analysis to assess differences between muscles was performed on the [increment]Ct. Both parametric (paired t tests) and non-parametric (Wilcoxon rank-sum test) tests were used to check the consistency of the results. The relative quantification between muscles was done using normalised equation M7 values (Livak and Schmittgen 2001).

Integration on the gene coexpression network

To infer coexpression patterns for the DE genes, the list of DE genes were input in the database of the bovine gene coexpression network (Reverter et al. 2006b). This gene coexpression network was constructed by integrating the data from 147 microarray hybridisations corresponding to nine seemingly independent experiments and spanning a total of 47 experimental conditions or treatments. Following normalisation via a multivariate mixed ANOVA model, significant gene-to-gene associations were identified using partial correlation coefficients coupled with an information theory approach. The numerical algorithm was later formalised (Reverter and Chan 2008). A sub-matrix of correlations between pairs of genes to depict coexpression patterns was obtained for the DE genes detected in this experiment that were present in the network.


Differentially expressed genes and validation by Taqman real-time RT-PCR

Model-based cluster analysis revealed that the model considering three clusters was the one that best fit the data. The central cluster contained those clones with the largest probability of being non-DE. Two FDR values, 5% and 10%, were imposed as cut-off points. For the FDR of 5%, 72 genes were identified as DE. After relaxing the FDR to 10%, 84 genes were deemed to be DE.

In order to validate the microarray results, differential expression of 11 DE genes was further analysed by TaqMan real-time RT- PCR. These genes were mostly a sample of the DE genes detected between the 5% FDR and the 10% FDR. Results were consistent both with the microarray experiment and between parametric and non-parametric tests, indicating that all genes showed significant differential expression between muscles. The p values were smaller or equal than 0.001 for all genes except for two, GREG1 and GPD1, which had p values smaller than 0.05 and 0.01, respectively. Supplementary Table S1 shows the log2 of the average n-fold differences in the muscle in which the corresponding gene appeared up-regulated (PM: ATF4, GPD1, MYH1, PDLIM7 and PFKM genes; FD: CA3, CREG1, CRYAB, CSRP3 and LDHB genes). AK1 did not show any level of expression in the FD muscle. Provided that all genes used in the validation experiment showed differential expression, the 84 genes found as DE considering a 10% FDR were further investigated.

Gene annotation and ontology

Putative gene functions could be assigned to 80 out of the 84 genes, 49 up-regulated in PM and 31 in FD. The remaining four genes, which are up-regulated in PM, matched to bovine sequences of unknown function. Table 1 provides a summary of the gene ontology (GO) annotation of the 80 DE genes. The gene symbol and name, the Gene ID, the number of DE array elements per gene and the normalised fold expression changes are shown in Supplementary Tables S2 and S3 for PM and FD, respectively. Most of the genes were represented by more than one clone on the microarray, ranging from 1 to 19, due to the redundancy of the non-normalised cDNA libraries used to build the array. The normalised fold changes ranged from 1.50 to 15.82, and for those genes represented by more than one clone, the minimum and maximum values are given.

The same GO functional groups appeared represented in both muscles, except for the extracellular matrix related group that was only represented in the FD muscle samples. The proportion of genes in each group varied between muscles. Thus, metabolic genes were three times more represented in PM (39%) than in FD (13%). The metabolic genes detected in PM were mainly involved in glucose and glycogen metabolic processes. On the other hand, genes related to muscle development and regeneration were three times more represented in FD (16%) than in PM (6%). The DE genes coding for contractile and constituent proteins of muscle, as well as those related to the different mechanisms involved in gene expression and protein modification, presented similar percentages in both muscles (24% and 31%).

Integration on the gene coexpression network

The list of 84 DE genes was input in the database of the bovine gene coexpression network (Reverter et al. 2006a) to infer coexpression patterns. As previously mentioned, this network has been constructed based on 47 experimental conditions and using the same microarray platform as the one used in the present study. The use of the same platform will provide a higher level of consistency in interpretation than would be achieved if different platforms were used.

Figure 1 shows the heat map of the coexpression correlation matrix for 61 out of the 84 DE genes. Thirty nine of them were up-regulated in PM (upper left corner of the correlation matrix) and 22 in FD (lower right corner). The remaining 23 DE genes did not appear in the coexpression network, as they were not differentially expressed in the experiments used to build it. Different groups of coexpression could be established according to the strong and positive correlation values (reddish colours) among genes in each muscle. The largest group was formed by 28 genes up-regulated in PM listed from ATP2A1 to MYOM2 in Fig. 1.

Fig. 1
Correlation matrix among the differentially expressed genes. The 39 genes located in the upper left corner of the correlation matrix corresponded to genes up-regulated in Psoas major. The 22 genes located in the lower right corner corresponded to genes ...


As in similar studies (Bai et al. 2003; Bernard et al. 2007; Campbell et al. 2001; Sudre et al. 2005), the microarray profiling of PM versus FD detected changes in gene expression that were mainly related to muscle contraction/structural constituents and to metabolic processes (54% of DE genes). Muscle contraction and constituents genes shown in Table 1 constituted one of the major groups of DE genes. MYH1, the gene coding for the IIX fibre type, was up-regulated in PM, the only muscle exhibiting the IIX fibres in this study. In addition to the myosin heavy chains, a profile for myosin light chains was observed. MYL1, MYL3 and MYLPF were up-regulated in PM, while MYL2 and MYL6 were up-regulated in FD. MYL1, MYL3 and MYLPF appeared in a gene coexpression group in Fig. 1 with MYH1, whereas MYL2 showed an independent pattern, which was consistent with its up-regulation in FD muscle. Thirteen genes coding for enzymes involved in the metabolism of sugars, Metabolic genes in Table 1, were up-regulated in PM, providing this muscle with the metabolic transcriptomic profile of a muscle with a predominance of fast fibre types. A group of six genes, from CKM to PYGM in Fig. 1, showed a strong positive coexpression pattern to genes from ACTN3 to MYH1, reflecting the physiological coordination required to store and sequester glycogen and phosphocreatine, the preferred substrates of the faster fibre types (Garrett and Kirkendall 2000). While glycolytic activity was significantly higher in PM than in FD, no differences were observed between the oxidative activities of PM and FD in previous experiments (Moreno-Sánchez et al. 2008), and the gene expression pattern agreed with those results. Besides the genes coding for glycolytic enzymes, the mitochondrial genes, ATP6, COX3 and CYTB, which have key roles in the oxidative phosphorylation pathway, were up-regulated in PM. Although these mitochondrial genes were not present in Fig. 1, CKMT2, a component of the mitochondria, was clustered in the coexpression group of the metabolic and structural genes and highly co-expressed with MYH1, reflecting the adaptations required for the oxidative characteristics of fast fibre types in PM.

Muscle specialisation is the result of the coordinated expression of contractile and metabolic proteins together with histological features that characterise the fibres (Acevedo and Rivero 2006; Quiroz-Rothe and Rivero 2004). The cellular presence of the various MyHC isoforms has been found to be coordinated with the metabolism (oxidative or glycolytic) and the morphological features in FD and PM (Moreno-Sánchez et al. 2008). As previously mentioned, expression of fast-glycolytic fibre related genes was positively correlated with expression of genes for glycolytic enzymes in PM (Fig. 1), which is consistent with the simultaneous up-regulation of metabolic and structural genes in PM. This indicates that the coordination established at the biochemical level also exists at the transcriptional level. However, the magnitude of the correlations between fibre attributes at both the biochemical (Moreno-Sánchez et al. 2008) and the transcriptomic levels (Fig. 1) suggested that properties of muscles may not be fully explained by their fibre type composition. Nevertheless, coordination of oxidative and contractile properties of FD could not be established because these two muscles did not exhibit oxidative differences and the difference in type I fibres between FD and PM (25%) was not supported at the transcriptional level. However, the gene encoding for the type I fibres (MYH7) consistently showed higher expression in FD than in PM (data not shown). Understanding these relationships will assist to ascertain the ability of muscles to react as a result of environmental factors such as exercise demand, changes in nutrition and disease challenge.

Organoleptic properties in meat may result from a combination of different biological processes (Bernard et al. 2007), and the relationship between meat quality and muscle attributes is not fully understood. The differential expression of some of the genes in Table 1 has been associated with different meat quality traits (Reverter et al. 2008). Specifically, the up-regulation of CRYAB and CSRP3 has been associated with low sensory scores of tenderness, flavour and juiciness (Bernard et al. 2007), which corresponded to the meat properties of FD in this comparison, in which CRYAB and CSRP3 were up-regulated. In addition, different heat shock proteins, including Hsp27 (LOC512251 in Table 1), have been described as down-regulated in distinct sensory traits (Bernard et al. 2007), favouring cytoskeleton actin disorganisation and degradation. Its up-regulation in FD might suggest that actin was mediating in the toughness of this muscle. Down-regulation of CRYAB, MUSTN1 and CSRP3 could also improve meat tenderness (Reverter et al. 2006b), agreeing with the fact that PM is more tender than FD (Díaz et al. 2006). Among other meat quality properties, differences in intramuscular fat content (IMF) have been described between FD and PM (Díaz et al. 2006), with PM showing higher IMF. However, two fat-related genes (FABP3 and CLU) appeared up-regulated in FD. In addition, up-regulation of CRYAB and FHL1 have been found to be associated with exercise endurance in mice (Wu et al. 2003); thus, their up-regulation in FD confirmed the endurance capability of this muscle. Although FD has an endurance-oxidative muscle profile, it requires further investigation as this profile was not supported by the DE of mitochondrial or oxidative metabolism related genes.

Changes in connectivity between genes across GO groups were found. Since the platform used in this experiment is the same one used to develop the network, it allowed a robust comparison with the large integrated analysis previously undertaken across a range of experimental conditions (Reverter et al. 2006a) and provided better consistency in interpretation than would we achieved if different platforms were used. Thus, four members of the troponin complex were DE, three fast forms (TNNC2, TNNI2 and TNNT3) in the PM and a slow form (TNNT1) in the FD (Table 1). Although all of them were involved in the striated muscle contraction and the transport of Ca2+, their coexpression patterns in Fig. 1 were different from the initially anticipated. The expression of TNNT1 was positively correlated to the rest of the troponins and genes up-regulated in PM, when it was expected to be independent or negatively correlated to them as they were up-regulated in different muscles. It represented a change in connectivity of these genes in the present comparison. Correlation patterns between troponins need further investigation as they differ in our experiment from the multiple conditions represented in the gene coexpression network (Fig. 1). In addition to these changes in connectivity (Reverter et al. 2006b) found between muscles, some others were also observed in the coexpression matrix, when the correlations among up-regulated genes within muscle are 0 or negative, as shown for FD in Fig. 1. Under the experimental conditions used to build the bovine network (Reverter et al. 2006a), the group genes from LDHB to MRCL3 in Fig. 1 were negatively correlated to the genes from KBTBD10 to UBA52. However, all of them were up-regulated in FD and, therefore, expected to show a positive correlation. It is noticeable that some of the genes up-regulated in FD had a different coexpression pattern under conditions such as nutritional restriction or different quality diets, with extracellular matrix (ECM) genes being down-regulated while CSRP3 and EEF1A1 were up-regulated (Reverter et al. 2004; Byrne et al. 2005; Lehnert et al. 2006). These genes appeared in many comparisons, and the changes in connectivity shown in our experiment illustrated how gene interactions can behave differentially under distinct experimental conditions.

Among the Transcription, RNA processing, protein synthesis and modification genes in Table 1, two transcription factors (TF), ATF4 and CREG1, were found up-regulated in PM and FD, respectively. Although recent studies in human did not find any association between TF mRNA levels and muscle type (Plomgaard et al. 2006), ATF4 correlated positively with all structural and metabolic genes up-regulated in PM (Fig. 1) and negatively or not correlated to ECM-related genes up-regulated in FD. Given that almost 20% of genes up-regulated in FD were related to ECM (Table 1), it is worth highlighting that the expression of these genes has been described to show a decreasing trend across muscle development (Hudson et al. 2009). Our results suggest a differential evolution of PM and FD. Differences in transcriptional regulation underpin much biological variation. TF would be the key in silencing genes related to ECM in PM and inducing them in FD muscle at adult ages. This point illustrated how the different specialisation and evolution of muscles can alter the established relationships between genes and open new hypothesis to a deeper investigation of the role of TF in muscle specialisation. Decoding causal transcriptional regulation remains challenging in muscle, as only a small number of well characterised TF is proposed to regulate development (Hudson et al. 2009). An additional change in connectivity affected a gene involved in translation in PM. RPS7, up-regulated in PM, showed a strong and negative correlation (blue colour) to structural and metabolic genes up-regulated in PM (Fig. 1). There were also some genes involved in translation up-regulated in FD (RPLP0 and EEF1A1), which, in agreement to our pattern of expression, showed a strong negative correlation to structural and metabolic genes up-regulated in PM. These results indicated that the ribosomal synthetic machinery was active in both muscles and that this up-regulation seemed to be coupled with the down-regulation of genes involved in fast-glycolytic features in the gene network conditions.

In conclusion, the present results confirm the existence of a muscle-dependent transcription pattern and the fact that the correlation among genes (which are not immutable) depends on the muscle type. Changes in connectivity were found in the different GO groups. These results also suggest the necessity of different muscle type data integration in order to perform a comprehensive network able to describe the transcriptional landscape of bovine skeletal muscle.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Table S1(46K, doc)

Primers and TaqMan MGB probes used in the validation by real-time RT-PCR, log2 of the average n-fold differences and statistical significance for every gene, except for the endogenous gene (ACTA1) and AK1, which was only expressed in FD (DOC 46 kb)

Supplementary Table S2(112K, doc)

Genes up-regulated in psoas major (DOC 111 kb)

Supplementary Table S3(71K, doc)

Genes up-regulated in flexor digitorum (DOC 70 kb)

Supplementary Fig. S1(2.5M, gif)

Experimental design for the microarray hybridisations. PM and FD refers to the Psoas major and Flexor digitorum muscles. Twenty slides were used to compare the PM and the FD muscles from ten animals. The sub-indexes indicate the animal sample, from 1 to 10. The arrows represent the hybridisations, and the red or green colours indicate the labelling of the sample (GIF 2604 kb)


The authors thank the Asociación de Avileña-Negra Ibérica and the Consejo Regulador de Carne de Ávila for their contribution in the sample collection. We are grateful to Dr. YongHong Wang for providing supports with the microarray methods and useful comments, to Dr. Jaime Cubero for skilled advice and the supply of the equipment for the real-time RT-PCR and to Dr. Sigrid Lehnert for lab support with the sequencing of some clones. We also thank Prof. Roger Cue (McGill University, Canada) for assistance in the English drafting of the manuscript. This work was partially supported by a grant from the Ministerio de Ciencia e Innovación of Spain (INIA-FEDER RTA2007-00081-00-00).

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.


  • Acevedo LM, Rivero JLL. New insights into skeletal muscle fibre types in the dog with particular focus towards hybrid myosin phenotypes. Cell Tissue Res. 2006;323:283–303. doi: 10.1007/s00441-005-0057-4. [PubMed] [Cross Ref]
  • Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–410. [PubMed]
  • Bai Q, et al. Development of a porcine skeletal muscle cdna microarray: Analysis of differential transcript expression in phenotypically distinct muscles. BMC Genomics. 2003;4:8. doi: 10.1186/1471-2164-4-8. [PMC free article] [PubMed] [Cross Ref]
  • Bernard C, et al. New indicators of beef sensory quality revealed by expression of specific genes. J Agric Food Chem. 2007;55:5229–5237. doi: 10.1021/jf063372l. [PubMed] [Cross Ref]
  • Brent AE, Tabin CJ. White meat or dark? Nat Genet. 2004;36:8–10. doi: 10.1038/ng0104-8. [PubMed] [Cross Ref]
  • Byrne KA, et al. Gene expression profiling of muscle tissue in Brahman steers during nutritional restriction. J Anim Sci. 2005;83:1–12. [PubMed]
  • Campbell WG, et al. Differential global gene expression in red and white skeletal muscle. Am J Physiol Cell Physiol. 2001;280:C763–C768. [PubMed]
  • Cassar-Malek I, et al. Muscle-specific metabolic, histochemical and biochemical responses to a nutritionally induced discontinuous growth path. Anim Sci. 2004;79:49–59.
  • Díaz C, Moreno-Sánchez N, Moreno A, Rueda J, Carabaño MJ (2006) Genetic basis of beef quality differences between muscles in beef cattle: Avileña negra-ibérica, a study case. In Proceedings of the XVI Congresso de Zootecnia, Castelo Branco (Portugal), 1–4 November, 21–26.
  • Díaz C, et al. Model selection in a global analysis of a microarray experiment. J Anim Sci. 2009;87:88–98. [PubMed]
  • Eden E, Navon R, Steinfeld I, Lipson D, Yakhini Z. Gorilla: a tool for discovery and visualization of enriched go terms in ranked gene lists. BMC Bioinform. 2009;10:48. doi: 10.1186/1471-2105-10-48. [PMC free article] [PubMed] [Cross Ref]
  • Efron B. Large-scale simultaneous hypothesis testing; the choice of a null hypothesis. J Am Stat Assoc. 2004;99:96–104. doi: 10.1198/016214504000000089. [Cross Ref]
  • Efron B, Tibshirani R. Empirical Bayes methods and false discovery rates for microarrays. Genet Epidemiol. 2002;23:70–86. doi: 10.1002/gepi.1124. [PubMed] [Cross Ref]
  • Flück M. Functional, structural and molecular plasticity of mammalian skeletal muscle in response to exercise stimuli. J Exp Biol. 2006;209:2239–2248. doi: 10.1242/jeb.02149. [PubMed] [Cross Ref]
  • Gardner TS, Bernardo D, Lorenz D, Collins JJ. Inferring genetic networks and identifying compound mode of action via expression profiling. Science. 2003;301:102–105. doi: 10.1126/science.1081900. [PubMed] [Cross Ref]
  • Garrett WE, Kirkendall DT (2000) Intramuscular energy stores and myoglobin. Exercise and Sport Science, pp. 78–80
  • Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using David bioinformatics resources. Nat Protoc. 2009;4:44–57. doi: 10.1038/nprot.2008.211. [PubMed] [Cross Ref]
  • Hudson NJ, Reverter A, Wang Y, Greenwood PL, Dalrymple BP. Inferring the transcriptional landscape of bovine skeletal muscle by integrating co-expression networks. PLoS ONE. 2009;4:e7249. doi: 10.1371/journal.pone.0007249. [PMC free article] [PubMed] [Cross Ref]
  • Kerr MK. Design considerations for efficient and effective microarray studies. Biometrics. 2003;59:822–828. doi: 10.1111/j.0006-341X.2003.00096.x. [PubMed] [Cross Ref]
  • Lehnert SA, Wang YH, Byrne KA. Development and application of a bovine cdna microarray for expression profiling of muscle and adipose tissue. Aust J Exp Agric. 2004;44:1127–1133. doi: 10.1071/EA03238. [Cross Ref]
  • Lehnert SA, Wang YH, Tan SH, Reverter A. Gene expression-based approaches to beef quality research. Aust J Exp Agric. 2006;46:165–172. doi: 10.1071/EA05226. [Cross Ref]
  • Lehnert SA, et al. Gene expression studies of developing bovine longissimus muscle from two different beef cattle breeds. BMC Dev Biol. 2007;7:95. doi: 10.1186/1471-213X-7-95. [PMC free article] [PubMed] [Cross Ref]
  • Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-δδct method. Methods. 2001;25:402–408. doi: 10.1006/meth.2001.1262. [PubMed] [Cross Ref]
  • Maccatrozzo L, Patruno M, Toniolo L, Reggiani C, Mascarello F. Myosin heavy chain 2b isoform is expressed in specialized eye muscles but not in trunk and limb muscles of cattle. Eur J Histochem. 2004;48:357–366. [PubMed]
  • McLachlan GJ, Bean RW, Ben-Tovim J, Zhu JX. Using mixture models to detect differentially expressed genes. Aust J Exp Agric. 2005;45:859–866. doi: 10.1071/EA05051. [Cross Ref]
  • Moreno-Sánchez N, Díaz C, Carabaño MJ, Rueda J, Rivero JL. A comprehensive characterisation of the fibre composition and properties of a limb (flexor digitorum superficialis, membri thoraci) and a trunk (psoas major) muscle in cattle. BMC Cell Biol. 2008;9:67. doi: 10.1186/1471-2121-9-67. [PMC free article] [PubMed] [Cross Ref]
  • Plomgaard P, et al. The mRNA expression profile of metabolic genes relative to MHC isoform pattern in human skeletal muscles. J Appl Physiol. 2006;101:817–825. doi: 10.1152/japplphysiol.00183.2006. [PubMed] [Cross Ref]
  • Quiroz-Rothe E, Rivero JLL. Coordinated expression of myosin heavy chains, metabolic enzymes, and morphological features of porcine skeletal muscle fiber types. Microsc Res Tech. 2004;65:43–61. doi: 10.1002/jemt.20090. [PubMed] [Cross Ref]
  • Reverter A, Chan EKF. Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks. Bioinformatics. 2008;24:2491–2497. doi: 10.1093/bioinformatics/btn482. [PubMed] [Cross Ref]
  • Reverter A, Byrne KA, Dalrymple BP (2003) Bayesmix: a software program for Bayesian analysis of mixture models with an application to model-based clustering of microarray gene expression data. In: XV Proc Assoc Adv Anim Breed Genet, Melbourne, VIC, Australia. pp. 90–93
  • Reverter A, et al. Joint analysis of multiple cdna microarray studies via multivariate mixed models applied to genetic improvement of beef cattle. J Anim Sci. 2004;82:3430–3439. [PubMed]
  • Reverter A, et al. Construction of gene interaction and regulatory networks in bovine skeletal muscle from expression data. Aust J Exp Agric. 2005;45:821–829. doi: 10.1071/EA05039. [Cross Ref]
  • Reverter A, et al. A gene co-expression network for bovine skeletal muscle inferred from microarray data. Physiol Genomics. 2006a;28:76–83. doi: 10.1152/physiolgenomics.00105.2006. [PubMed] [Cross Ref]
  • Reverter A, Ingham A, Lehnert SA, Tan SH, Wang YH, Ratnakumar A, Dalrymple B. Simultaneous identification of differential gene expression and connectivity in inflammation, adipogenesis and cancer. Bioinformatics. 2006b;22:2396–2404. doi: 10.1093/bioinformatics/btl392. [PubMed] [Cross Ref]
  • Reverter A, et al. Dissection of beef quality phenotypes using a myogenin network-anchored systems biology approach. Aus J Exp Agric. 2008;48:1053–1061. doi: 10.1071/EA08052. [Cross Ref]
  • Segal E, et al. Module networks: identifying regulatory modules and their condition-specific regulators from gene expression data. Nat Genet. 2003;34:166–176. doi: 10.1038/ng1165. [PubMed] [Cross Ref]
  • Stuart JM, Segal E, Koller D, Kim SK. A gene coexpression network for global discovery of conserved genetic modules. Science. 2003;302:249–255. doi: 10.1126/science.1087447. [PubMed] [Cross Ref]
  • Sudre K, et al. Transcriptome analysis of two bovine muscles during ontogenesis. J Biochem (Tokyo) 2003;133:745–756. [PubMed]
  • Sudre K, et al. Biochemical and transcriptomic analyses of two bovine skeletal muscles in charolais bulls divergently selected for muscle growth. Meat Sci. 2005;70:267–277. doi: 10.1016/j.meatsci.2005.01.012. [PubMed] [Cross Ref]
  • Tan SH, Reverter A, Wang YH, Byrne KA, McWilliam S, Lehnert SA. Gene expression profiling of bovine in vitro adipogenesis using a cdna microarray. Funct Integr Genomics. 2006;6:235–249. doi: 10.1007/s10142-005-0016-x. [PubMed] [Cross Ref]
  • Tanabe R, Muroya S, Chikuni K. Sequencing of the 2a, 2x, and slow isoforms of the bovine myosin heavy chain and the different expression among muscles. Mamm Genome. 1998;9:1056–1058. doi: 10.1007/s003359900924. [PubMed] [Cross Ref]
  • Tanner CJ, et al. Muscle fiber type is associated with obesity and weight loss. Am J Physiol Endocrinol Metab. 2002;282:E1191–E1196. [PubMed]
  • Thompson J. Managing meat tenderness. Meat Sci. 2002;62:295–308. doi: 10.1016/S0309-1740(02)00126-2. [PubMed] [Cross Ref]
  • Wang YH, et al. Transcriptional profiling of skeletal muscle tissue from two breeds of cattle. Mamm Genome. 2005;16:201–210. doi: 10.1007/s00335-004-2419-8. [PubMed] [Cross Ref]
  • Wang YH, et al. Gene expression patterns during intramuscular fat development in cattle. J Anim Sci. 2009;87:119–130. doi: 10.2527/jas.2008-1082. [PubMed] [Cross Ref]
  • Wu H, Gallardo T, Olson EN, Williams RS, Shohet RV. Transcriptional analysis of mouse skeletal myofiber diversity and adaptation to endurance exercise. J Muscle Res Cell Motil. 2003;24:587–592. doi: 10.1023/B:JURE.0000009968.60331.86. [PubMed] [Cross Ref]

Articles from Springer Open Choice are provided here courtesy of Springer