Abstract: Gene expression in different growth stages showed differences. In order to establish the correlation between gene expression mechanisms and beef traits formation which controlled by polygenes, microarray analysis was performed in Longissimus dorsi muscle tissue between 1 and 24 months Chinese red steppes. About 1,282 (5.6%) probes showed significant differences at the 2 growth stages and 126 genes showing strong correlation with beef traits formation were gained by the GO analysis. With the KEGG analysis, 63 pathways were found to be related to beef traits formation which involved 73 genes. About 28 genes were found in a single pathway while 35 genes were found in 2-16 pathways, respectively. The panel of transcripts and gene pathway analysis in different growth stages may be helpful for the study on beef traits formation and the gene expression profile construction in Longissimus dorsi muscle tissues would also make a model for screening candidate genes which have genetic effect on meat quality in bovine.

INTRODUCTION

Beef traits controlled by polygenes are affected by many factors such as genetic factors, nutrition, age, slaughtering and processing technology but the genetic factors are the most influential. In-depth studies of the molecular mechanisms of beef traits formation are expected to develop. Kee et al. (2008) analyzed the correlation between beef transcriptome and 4 different beef traits: Shear Force (SF), Water Holding Capacity (WHC), Cooking Loss (CL) and Loin Eye Area (LEA) in Korean native cattle. The results showed that SF was mainly related to energy metabolism, LEA to rRNA processing, WHC to protein metabolism. Thus energy, protein metabolism and meat quality establishes a causal relationship after the aging process. In Ponsuksili et al. (2008), the resarcher showed that in 1,279 transcripts with trait correlated expression to WHC, 897 expressions QTL (eQTL) with 104 eQTL coinciding with QTL regions for WHC, 96 transcripts had trans acting and 8 transcripts had cis acting regulation. These results showed that the differences of traits may eventually be linked to genes expressions. Thus, we may predict the different traits by analyzing gene expression mechanisms. Te Pas et al. (2007) showed that 88 pathways of porcine myogenesis are involved by Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) analysis. In their research, most genes were found in a single pathway, even some were found in up to seven pathways. Combining the pathways and the microarray information, about 21 pathways showed sufficient information content for further analysis. These pathways were related to the regulation of several steps in myogenesis and energy metabolism. The experiment indicates that combining microarray results and biological pathway information will give insight into biological trait formation.

The present study try to relate gene expression changes to beef traits formation at two stages (1 and 24 months of age) by the GO and the KEGG analysis. Knowing these relationships may provide a better understanding of gene expression in beef traits controlled by polygenes. It was necessary to develop method to select candidate genes of beef traits by searching Gene Ontology (GO) and biochemical pathways databases (KEGG).

MATERIALS AND METHODS

Animals and muscle samples: A total of 18 bulls of the same breed (Chinese red steppe) were included in the study; 9 were 1 month old while remaining 9 were 24 months old and all provided by Jilin Academy of Agricultural Sciences. They were maintained in standard conditions and fed with standard diets. The bulls were humanly killed at the slaughter house of the academy and fresh longissimus dorsi muscle tissues were obtained during the process which were immediately frozen in liquid nitrogen and stored at -80°C for microarray analysis.

Total RNA isolation and validation: Total RNA from muscle samples was extracted using Trizol reagent (Invitrogen, USA) and purified using QIAGEN RNeasy Total RNA Isolation kit (QIAGEN, USA). The RNA quantity was assessed by measure the optical density at 260/280 nm using Nanodrop ND-1000 (Nano-Drop Technology, Wilmington, DE) and the integrity of total RNA samples were assessed qualitatively on an agilent bio analyzer by inspection of 18 and 28S rRNA bands after gel electrophoresis.

Microarray target preparations and hybridization: The microarray analysis (probe labeling, hybridization and scanning) was performed using Affymetrix GeneChip Bovine Genome Array (Affymetrix, Santa Clara, CA, USA) following the manufacturer’s instruction in the CapitalBio Corporation, Beijing, China. Total RNA (10 μg) or mRNA (0.2 μg) was first reverse transcribed using a T7-Oligo (dT) promoter primer in the first-strand cDNA synthesis reaction. Following RNase H-mediated second-strand cDNA synthesis, the double-stranded cDNA was purified and served as a template in the subsequent In vitro Transcription (IVT) reaction. The IVT reaction was carried out in the presence of T7 RNA Polymerase and a biotinylated nucleotide analog/ribonucleotide mix for complementary RNA (cRNA) amplification and biotin labeling. The biotinylated cRNA targets were fragmented in 1 x fragmentation buffer solution provided with the genechip sample cleanup module (Affymetrix) at 94°C for 35 min. A total of 10 μg of fragmented biotin-labeled cRNA per replicate in hybridization mixture then was hybridized to Bovine Genome Array from Affymetrix GeneChips and incubated overnight at 45°C in Affymetrix GeneChip Hybridization Oven 640, all according to the manufacturer’s instructions. The mixture was removed 16 h after hybridization in several cycles; the chips were washed with nonstringent buffer and stained with streptavidin-phycocrythrin antibody solution (Affymetrix) on an automated Affymetrix GeneChip Fluidic Station 450 station. The data were collected by using affymetrix genechip scanner 3000. Microarray images quantified using Affymetrix GeneChip Operating Software (GCOS1.4).

Quantitative real-time RT-PCR: To validate the results of the microarrays, 6 genes (IGF2 (down-regulated), PSPH (down-regulated), IGFBP2 (no significant change), IGFBP5 (no significant change), GPAM (up-regulated), CAPN3 (up-regulated)) showing differential expression were selected and analyzed by real-time RT-PCR. Genes were chosen in different pathways. Primers were designed using Primer Premier5 software (Table 1). Relative mRNA expression levels were determined by using the cDNA for quantitative real-time PCR amplification with SYBR Green I kit (Biotechs, China). Specificity of all individual amplification reactions was confirmed by melting curve analysis. Real-time expression values were calculated by using the relative standard curve method. Standard curves were generated for each mRNA by using 10-fold serial dilutions. The reactions were denatured at 95°C for 5 min followed by 45 cycles at 95°C for 10 sec and 60°C for 30 sec in the ABI PRISM 7500 Sequence detection system (Applied biosystems, USA). Data are reported as values normalized to the housekeeping gene GAPDH. One way ANOVA analysis was performed using SPSS (13.0) software. The significance of differences between the means for each gene was evaluated using LSD or Tamhane’s T2 test (for significance levels of p<0.05). All the values are presented as means±SEM.

Table 1:	Primer pairs used for real-time PCR

Microarray data analysis: The 6 chips rectified were made by normalization processing using dChip (Affymetrix) software (Li and Wong, 2001a, b). False discovery rates for the genes were calculated by using t-test and p-values. Fold changes were calculated based on unadjusted data’s means.

Genes were annotated with NetAffx (Liu et al., 2003) and differentially expressed genes were identified using Significance Analysis of Microarrays (SAM). Using the criterion of cutoff limitation as a fold change2 or 0.5 and q-value 5%, differential expression genes were screened and clustered with the Cluster 3.0 software. The selected genes were further analyzed in the context of Gene Ontology (GO) biological process and KEGG biological pathway using the molecule annotation system 2.0 (MAS 2.0, http://bioinfo.capitalbio.com/mas3/) software (CapitalBio, Beijing, China).

RESULTS AND DISUSSION

Microarray analysis: The Affymetrix GeneChip Bovine Genome Array which contains 24,128 probe sets representing >23,000 transcripts and 19,000 UniGene clusters was used to profile gene expression in the longissimus dorsi muscle of 1 and 24 months old Chinese red steppes.

The microarray hybridization results showed that 22,810 probe sets were used to detect which accounts 94.5% of the total number of probes and 16,321 (71.6% of the total number) probe sets were detected. At the cutoff criteria of the signed fold change 2 or 0.5 and qQ5%, a total of 1,282 transcripts (5.6% of total probe sets on the array) were differentially expressed at the two stages, 959 (4.2%) of which were up-regulated while 323 (1.4%) were down-regulated (GEO number: GSE21782). About 1,008 differential gene annotations obtained using the molecule annotation system 2.0 (CapitalBio, Beijing, China) among of which 4 gene transcripts are unknown proteins in the functions (0.4%) while 274 of which show no homology to any known sequence.

BLAST analysis revealed that 1,001probe genes (745 up-regulated and 256 down-regulated) share significant similarity in amino acid sequence with other functional genes accounting for 78.1% of the total number of probe sets. However, 281 (21.9%) probe genes (214 up-regulated and 67 down regulated) did not show similarity with other functional genes (June 2009).

Quantitative real-time RT-PCR: To confirm the microarray data, 6 genes (2 down-regulated genes (IGF2 and PSPH), 2 up-regulated genes (GPAM and CAPN3) and 2 genes with no significant change (IGFBP2 and IGFBP5)) with differentially expressed were selected for quantitative real-time PCR. Absolute quantification method was chosen for each of the 6 genes in the experiment (Tichopad et al., 2003). Linearized plasmid DNA was used to plot a standard curve to quantify absolute gene expression changes and the R² for each standard curves was >0.99. The values were presented as means±SEM (Table 2). SPSS analysis showed that the expression of IGF2 (Sig., 0.008) and PSPH (Sig., 0.032) are down-regulated, GPAM (Sig., 6.919e-008) and CAPN3 (Sig., 0.004) are up-regulated, IGFBP2 (Sig., 0.808) and IGFBP5 (Sig., 0.748) with no significant change (Fig. 1). These results were exactly consistent with the microarray experiment.

Gene Ontology (GO) analysis: The selected genes with significant expression difference at the 2 stages were analyzed in the context of GO biological process. The GO analysis dynamically assigns different genes to different GO biological process categories and provides a z score for each category.

The z score represents a statistical measurement of relative gene expression dynamic in a given GO term. This approach revealed very interesting phenomena. When we choose GO class of level 0, gene expression fall in 4 Gene ontology categories including molecular functions (464), biological processes (747) ontology categories and cellular components (201) but when it comes to the GO class of level 1,282 categories were obtained (Fig. 2).

Analysis of the genes may affect beef traits: About 585 genes in 1,282 genes of significant expression difference were provided with gene annotations in relation to cellular component. About 126 genes related to cellular component and metabolism out of 457 genes with definite biological process showed correlation with meat quality trails.

Based on the functional analysis of 126 genes showing correlation with beef traits, we identified a few common biological metabolism categories including glucose (16 genes, 10 up-regulated and 6 down-regulated), lipid (19 genes, 17 up-regulated and 2 down-regulated), protein and amino acid (40 genes, 32 up-regulated and 8 down-regulated), DNA and RNA (4 genes, 1 up-regulated and 3 down-regulated), energy (9 genes, 6 up-regulated and 3 down-regulated) and others (38 genes, 30 up-regulated and 8 down-regulated) (Table 3).

Table 2:	Number of copies of mRNA in different genes (copies μL-1)
(p<0.05)

KEGG analysis showed 63 pathways related to metabolism. About 73 genes expressed with significant difference. We also find 28 genes in a single pathway, 35 genes in 2-16 pathways (Table 4).

Meat quality traits were quantitative traits controlled by polygenes. Due to the special characteristic, many of genes could not be clearly established links with traits, although they played an important role. Some of them regulating the expression of other genes could not simply be correlated with in the way of measuring traits and gene expression. Therefore, measurement seemed to be pointless in this case. Thus, new methods were expected to develop.

Fig. 1:	The gene expression of IGF2, PSPH, IGFBP2, IGFBP5, GPAM, CAPN3 of 1 and 24 months Chinese red steppes longissimus dorsi muscle in quantitative real-time RT-PCR

Fig. 2:	The results of Gene Ontology (GO) analysis

Table 3:	Genes of significant differential expression related to beef traits in microarray analysis

Table 4:	The KEGG analysis on genes of significant differential expression related to beef traits in microarray analysis
*N = The number of genes with significant differential expression

In this study, microarray analysis was carried out to profile gene expression in Chinese red steppes Longissimus dorsi muscle at two different stages: 1 and 24 months old (the experiments have been achieved recognition of GEO and successfully obtained access number. Fortunately, 1,282 (5.6%) probes showed significant differences between the two growth stages in bovine and 126 genes among of which showing strong correlation with beef traits’ formation were gained by GO and KEGG analysis. These genes related to various metabolisms may be the reason which accounts for the differences of beef traits between the two stages. If we put the cutoff limitation scale down to a fold change 1.5 or 0.67 and q-value 5%, 2,923 (12.8%) probe sets were selected but in the case of the majority, the criterion of cutoff limitation may be not strict enough.

We also analyzed several genes which thought to be important and influential beef traits such as those encode Ca²⁺/calmodulin-dependent protein Kinases (CaMKs), myogenic regulatory factors (MRFs), Fatty Acid-binding Proteins (FABPs) and Peroxisome Prolifemtor-activated Receptors (PPARs). Calpains are a family of calcium activated cystein proteases that are widely distributed in cell cytoplasm and sub-organelles.

They go along with their endogenous inhibitor, calpastatin. In calpain families, calpain l-10 are in relation to tenderization of meat (Kar et al., 2010). In this study, calpain 1 (fold change = 0.97, q-value = 46.28%) has no significant changes. This indicates that calpain 1 may be unresponsible to the differences of beef traits in the experiment. Meanwhile, the experiment shows that Calpastatin (fold change = 2.13, q-value = 0.07%) was significantly up-regulated which indicates that the differences may be due to Calpastatin. The experiment also shows that calpain 3 was significantly up-regulated but in literature there are controversies on whether calpain3 has any effect on tenderization of meat. Correlation analysis between the relative rate of tenderization and mRNA expression of calpains revealed a strong relationship with calpain 3 in both species, bovine and ovine [r (bovine) = 0.522, r (ovine) = 0.706] (Ilian et al., 2001). In contrast, Parr et al. (1999) did not find a link between calpain 3 and postmortem tenderization in the porcine Longissimus Thoracis et lumborum (LT).

The discrepancy between them may be due to the methodology or animal species. In addition in this research, the calpain 10 also show no significant changes. Evidence is mounting that the Ca²⁺/calmodulin-dependent Kinases (CaMK) II and IV play important role in regulating oxidative enzyme expression, mitochondrial biogenesis and expression of fibre-type specific myofibrillar proteins.

Examination of Ca²⁺-dependent and Ca²⁺-independent CaMK II activity in chicken anterior latissimus dorsi muscle sample after 7 days of stretch overload and 2 weeks of voluntary wheel running have shown that stretch overload and wheel running result in a 122% increase in muscle protein content and are associated with 47% increase in Ca²⁺-independent CaMKII activity (Fluck et al., 2000).

Results indicate that CaMK IIg not CaMK IIa or -b is up-regulated in aged and denervated soleus muscle and CaMK IV is absent in skeletal but not cardiac muscle (Chin, 2004). SiDe milarly this study showed that CaMK Iid (fold change = 2.97, q-value = 0) and CaMK IIg (fold change = 3.06, q-value = 0) have significant changes in contrast with CaMK IIa, CaMK IIb or other members of CaMKs. Thus, the differences in beef traits may be related to the regulations of CaMK IId or CaMK IIg.

The MyoD gene family members encoding basic helixloop-helix (bHLH) transcription factors (MyoD, myogenin, Myf5 and MRF4) contain one or two transactivation domains (at the N- and C-terminus). They are able to make nonmuscle cells to have muscle phenotype in culture (Wei and Paterson, 2001). It is reported that MRFs have relation to tenderness, Loin Eye Area (LEA) and other beef traits (Kim et al., 2009; Liu et al., 2008; Zhu and Li, 2005), however the experiment showed that MyoD, MyoG, Myf6 and Myf5 have no significant changes. For MyoD (fold change = 1.33, q-value = 5.31%), although the q-value (5.31%) was in accordance with the criterion of cutoff limitation, the fold change (1.33) was not within the scope of screening.

So, it showed that the result of gene expression comparison has no significant difference even if the relaxation of screening criteria to fold change is relaxed to R1.5 or Q0.67 and q-value Q5%. So, the MyoD gene expression showed that differences in beef traits may be unrelated to regulation of MRFs.

Intramuscular fat content (IMF) was correlated with meat quality trails and showed significant effects on three aspects of meat: tenderness, flavor and juicy (De Vol et al., 1988). Fatty Acid-binding Proteins (FABPs) are small intracellular proteins involved in fatty acid transport from the plasma membrane to the sites of β oxidation and/or triacylglycerol or phospholipid synthesis (Veerkamp and Maatman, 1995). The experiment showed that FABPs have no significant changes even the relaxation of screening criteria to fold change is relaxed as above. So, the IMF expression results showed that the differences in beef traits may be unrelated to the regulation of FABPs.

The Peroxisome Proliferator-activated Receptors (PPARs) are members of the nuclear receptor superfamily and play an important regulatory role in lipid metabolism (Latruffe and Vamecq, 1997). The PPAR gamma (PPARG) subtype is a key transcription factor that controls adipocyte differentiation and fat deposition in mammals but in this experiment PPARG (fold change = 0.96, q-value = 46.28%) as well as other members of this family have no significant changes. The PPAR expression results showed that the differences in beef traits may be unrelated to the regulation of PPARs.

CONCLUSION

In this study, we have shown several CaMKs encoded genes play an important role in beef traits using Affymetrix GeneChip Bovine Genome Array while genes encoding MRFs, FABPs and PPARs do not have influential impact on beef traits’formation. We also found, it was an effective way to exclude invalid genetic pathways when the subject (beef traits in the model) was controlled by various genes. It is no doubt that we need to improve the efficiency of searching for candidate genes. At the same time, we also see the interaction of gene pathways and analyze the mechanisms of gene expression may provide some new insight in regulatory network characters of quantitative traits controlled by polygenes and the result will provided in due course.

ACKNOWLEDGEMENTS

We thank Melvin Korkor (Jilin University, P.R. China) and Naizhong Cui (Dalian University of Technology, P.R. China) for critical reading the manuscript, Yonghe Ma (Jilin University, P.R. China) for technical help and Supertech (Enuo Technology, P.R. China) for revision. This research was supported by the 863 hightech project (2008AA101010) and the National Scientific and Technological Supporting Project (2007BAD55B01).