Transcriptional profiling identifies differentially expressed genes in developing turkey skeletal muscle

Animals and Tissue Collection

Animals used in this study were domestic turkeys maintained at the Poultry Research Center at The Ohio Agricultural Research and Development Center located in Wooster, OH under AGACUC protocol #04-AG007. These turkeys were from two genetic lines: RBC2, a randomly bred control line representative of a 1967 commercial turkey [10, 19], and F, a subline genetically selected from the RBC2 line for increased 16wk body weight [11]. The birds were euthanized immediately prior to tissue harvest at one of three developmental stages: 18 d embryo (18de; hyperplasia), 1 d post hatch (1d; shift from myoblast-mediated growth to satellite cell modulated growth by hypertrophy), and 16 wk old birds (16wk; approximate age of commercial slaughter). Pectoralis major muscle was excised [17], snap frozen in liquid nitrogen, and stored at -80°C until shipment on dry ice to Michigan State University (MSU), where samples were stored at -80°C until RNA extraction.

RNA Extraction

Total RNA was isolated from skeletal muscle tissue using TRIReagent (Molecular Research Center, Inc., Cincinnati, OH) as per manufacturer instructions. All samples were treated with DNA-free™ DNase (Ambion, Inc., Austin, TX) to remove any contaminating genomic DNA. Sample integrity was confirmed using an Agilent 2100 Bioanalyzer (Santa Clara, CA) that determined RNA Integrity Number (RIN). All samples used for microarray analysis had an RIN > 8.0.

Microarray Experimental Design

The TSKLMO array was designed and created based on sequences from cDNA libraries constructed at crucial stages in muscle development [18]. Description of this platform can be found on the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database [GEO: GPL9788]. Experiment 1 (Additional File 1, Figure S1a) was designed to directly compare the three developmental stages, 18de vs. 1d or 1d vs. 16wk, for each of the genetic lines. Experiment 2 (Additional File 1, Figure S1b) was designed to directly compare the two genetic lines, F vs. RBC2, at each of the developmental stages. Each Experiment contained 10 arrays per comparison; therefore, Experiment 1 contained 40 arrays, and Experiment 2 contained 30 arrays. Birds were randomly assigned to array, and hybridizations were performed in random order.

RNA Amplification and Microarray Hybridization

Total RNA for microarray analysis was amplified and dye-coupled in preparation for hybridization using the Amino Allyl MessageAmp™ II aRNA Amplification Kit (Ambion, Inc.) per manufacturer instructions. All microarray preparation and hybridization procedures were performed as previously described [17]. Briefly, all microarrays were UV-crosslinked and pre-hybridized prior to hybridization. Fragmented, Cy3-coupled aRNA was mixed with its Cy5-coupled partner, and the mixtures were hybridized to oligonucleotide microarrays for 18 h in a GeneTac Hybridization Station (Genomic Solutions, Ann Arbor, MI). Arrays were then rinsed, dried by centrifugation, and immediately scanned with a GenePix 4000B scanner (Molecular Devices, Sunnyvale, CA). Image analysis was performed using GenePix Pro 6.0 software, and spot intensities were exported as GPR files for statistical analysis. Fluorescence intensity data were not background-corrected in accordance with recent recommendations [20, 21] but were normalized for dye intensity bias using the LOESS procedure of the normalizeWithinArray function of the Bioconductor R software LIMMA [22].

Statistical Analysis

The two sets of experiments were analyzed separately. For the first experiment, in which 1d was directly compared against the other two stages within arrays and within lines, all microarray data were expressed as log₂ ratios of each stage with 1d for each array. These log₂ ratios were analyzed using a linear mixed model based on an overall intercept, the effect of dye assignment for samples deriving from stage 1d, effects of the other two stages specified as numerator terms in the log₂ ratios, effect of line, and the random effects of the animal assigned to stage 1d and its interaction with the other two stages to serve as the experimental error term for stage as each oligonucleotide was spotted twice on an array. For this experiment, the comparisons of interest included all stages against each other within each of the two lines, as well as overall differences between the three stages.

For the second experiment, in which lines were compared directly against each other within arrays and within stages, all data were expressed as log₂ ratios of F over RBC2 for each array. These log₂ ratios were analyzed using a linear mixed model based on an overall intercept, the dye assignment for samples deriving from line RBC2, stage effects, the effect of line, and the random effect of the animal derived from the RBC2 line. The random effect could have been also arbitrarily chosen as the animal derived from the F line, but, in either case, this specification ensured that replication was based on the number of animals and not on the number of spots per gene per array. In this experiment, the comparisons of interest were overall line differences and line differences within each stage. SAS PROC MIXED was used for both analyses given the linear mixed model specifications for both experiments [23]. Statistical significance was based on the estimated false discovery rates (FDR) [24].

Raw Cy5 intensities, raw Cy3 intensities, LOESS-normalized log₂ Cy5:Cy3 ratios, and LOESS-normalized log₂ average intensities (A), as well as original GPR files for all arrays, were submitted to the National Center for Biotechnology Information (NCBI)'s Gene Expression Omnibus (GEO) [GSE19585]. Control probes were also evaluated for both Experiments 1 and 2. Results examining the fluorescence intensities of probes designed to contain an increasing number of mismatched bases as well as scrambled sequence oligos that served as negative controls indicated that hybridization was specific in these experiments. These data have been previously presented [17].

Gene annotation and functional analysis

After statistical analysis, genes were sorted by FDR, and the top 200 genes differentially expressed by overall effect of development (Experiment 1) and all genes differentially expressed by overall effect of genetic line (Experiment 2) were investigated by subjecting oligonucleotide sequences to an NCBI BLAST search and choosing best matches for accession number and putative annotation. In addition, the top 30 up-regulated and down-regulated genes in Experiment 1 and all differentially expressed genes in Experiment 2 in each comparison: F 18de:1d, F 1d:16wk, RBC2 18de:1d, and RBC2 1d:16wk were annotated in the same manner. These lists were then probed for patterns in expression using the Ingenuity Pathways Analysis software (IPA; Redwood City, CA).

Quantitative real-time PCR confirmation

Differential expression of skeletal muscle genes

Functional and pathway analysis

qPCR confirmation of differentially expressed genes

Genes in the categories of ECM function, cell death, and skeletal muscle development and function, along with several genes with miscellaneous functions, were selected for qPCR confirmation of microarray results. In addition, genes common to both Experiments 1 and 2 were selected for confirmation.

Figure 1 shows results for versican (1A), glypican-1 (1B), syndecan-4 (1C), spondin 2 (1D), fibroblast growth factor 2 (FGF2; 1E), and collagen VI α1 (COL6A1; 1F), which are all involved in ECM regulation. The COL6A1 gene was also common between Experiments 1 and 2. Versican and glypican-1 appeared to be expressed predominantly at the 18de stage when compared to later stages in development. Interestingly, spondin 2 was expressed almost exclusively at the 1d stage, during muscle hypertrophy. Several differences between the RBC2 and F lines were observed by qPCR analysis that were not detected with the significance threshold used for the microarray analysis. Most notably, FGF2 expression at the 18de stage was twice as high in the RBC2 line compared to that of the F line and nearly 3 times greater at the 1d stage in the RBC2 line compared to that of the F line.

Figure 2 shows results for genes involved in cell death/apoptosis as determined by IPA: caspase-3 (2A), growth arrest-specific 2 (GAS2; 2B), death-associated protein (DAP; 2C), β-parvin (2D), platelet-derived growth factor α (PDGFα; 2E), nestin (2F), matrix γ-carboxyglutamic acid (GLA) protein (MGP; 2G), and gelsolin (2H). The gene β-parvin was also common to both Experiments 1 and 2. Caspase-3, GAS2, DAP, and β-parvin all appeared to be down-regulated as muscle develops. Nestin was highly expressed at the 18de stage with greater than 100 times the expression of subsequent developmental stages. The MGP gene, however, was expressed principally at the 1d stage at 80.28 and 55.38 times the 18de expression in RBC2 and F lines, respectively. Gelsolin expression appeared to increase with muscle maturation and peak at the 16wk time point.

Genes important to skeletal muscle growth and development and calcium signaling are shown in Figure 3. These are myostatin (3A), Activin IIB receptor (ActIIBR; 3B), calmodulin (3C), calreticulin (3D), and isotypes of troponin: TNNC1 (3E), TNNI1 (3F), TNNT2 (3G), and TNNT3 (3H). Myostatin expression dipped at 1d posthatch, when hypertrophy occurs in muscle; interestingly, its expression was higher in the RBC2 line at 18de but was higher in the F line by 16wk of development. In conjunction, the myostatin receptor ActIIBR was highly down-regulated at 16wk (18de expression was 18.69 and 10.95 times greater in RBC2 and F lines, respectively), when myostatin expression is high. Relative calmodulin and calreticulin expression tended to respond in opposite directions with calmodulin expression highest at the 16wk stage and calreticulin highest at 18de and decreasing as muscle matures. Of the troponins, TNNI1 and TNNT2 exhibited similar profiles with predominant expression at 18de that decreased substantially at the other stages. Interestingly, TNNT3 expression appeared to be regulated in the opposite direction of TNNT2.

The expression profiles of several other genes that were common to both Experiments 1 and 2 but did not easily fall into a functional category, along with one gene with unknown function, are shown in Figure 4. These genes included: BTBD11 (4A), cMet (4B), plastin3 (4C), 5'-nucleotidase, cytosolic III (NT5C3; 4D), and nebulin-related anchoring protein (NRAP; 4E). In addition, a gene with no known annotation was chosen from the top 40 genes differentially expressed by overall effect of development (FDR < 1.0 × 10^-13) and was subsequently coded Unknown (4F). Expression of these genes was variable although line differences were observed at 1d for all genes. Expression of NRAP was dramatically up-regulated at 1d and 16wk when compared to 18de. Expression profiles for these 28 genes confirmed that microarray results were repeatable.

Extracellular Matrix Genes

Previous studies have demonstrated that the ECM comprising the connective tissue surrounding muscle can interact with growth factors, regulate cellular signal transduction pathways, and affect growing and developing muscle fibers [34–36]. Genes that fell into numerous functions or canonical pathways as well as some genes that were not included in IPA networks were selected for this group. The ECM is composed mainly of collagens and non-fibrous glycoproteins like proteoglycans. The expression of three proteoglycans: versican, glypican-1, and syndecan-4, was quantified in this study. Membrane-associated proteoglycans in the ECM are known to regulate numerous growth factors, the major constituents that regulate activation, proliferation, and differentiation of satellite cells [8]. A predominant example is FGF2, a potent stimulator of muscle fiber proliferation and a strong inhibitor of differentiation [35, 37]. This function would explain the high FGF2 expression at 18de and 1d in the current study, as these are periods of myoblast or satellite cell proliferation. The interaction between FGF2 and the heparan sulfate proteoglycans glypican-1 and syndecan-4 is required for high-affinity binding of FGF2 to its cellular receptor [35]. Differences in expression of syndecan-4 have been observed in vitro between developing RBC2 and F line satellite cells [15], and RBC2 satellite cell expression of syndecan-4 as well as glypican-1 were altered with the addition of FGF2 [38]. Satellite cells from syndecan-4 knockout mice display defective patterns of activation, proliferation, and differentiation [39], indicating a crucial role in muscle development that is also supported for turkeys in the current study. In addition, syndecan-4 and FGF2 appear to have similar expression profiles, especially at 1d. Versican is a chondroitin sulfate proteoglycan that was expressed predominantly at the 18de stage in this study. Fernandez et al. (1991) hypothesized that this protein's high embryonic expression at the time of fiber formation may be involved in fiber spacing, establishing the morphological structure of the muscle in chicks [40]. The space between muscle fibers may play a critical role in the water-holding capacity of mature muscle and may therefore be associated with quality defects like PSE meat. It is therefore possible that characteristics that contribute to meat quality are decided very early in myogenesis.

Dominant negative mutations in the COL6A1 gene, the only collagen that was investigated in the current study, are known to lead to Bethlem myopathy [41] and Ulrich congenital muscle dystrophy [42]; COL6A1 knockout mice exhibit sarcolemmal disorganization [43]. Taken together, it is clear that expression of genes whose proteins are located in the ECM of skeletal muscle likely contribute significantly to myogenesis in RBC2 and F turkeys.

Cell Death/Apoptosis Genes

Cell death and apoptosis are necessary events in muscle differentiation and maintenance. In addition, other reports have shown similarities between mechanisms of apoptosis and myoblast differentiation, suggesting that similar pathways and effectors may be in use [44, 45]. Therefore, it was logical to further evaluate the expression of several differentially expressed genes associated with these functions. For example, caspase-3, an "executioner" cysteine protease in apoptosis [46], is necessary for proper differentiation of C2C12 myoblasts [47]. However, several reports have shown that this role is separate from apoptosis [44, 45] although it employs similar mechanisms, such as DNA strand breakage [48].

A component of the microfilament, GAS2, is a substrate for caspase-3 and is involved in microfilament reorganization, possibly an early event in apoptosis [49]. The protein is expressed in embryonic mouse limbs and involved in apoptosis during interdigital development; GAS2 may also be involved in myogenesis, possibly in the fusion of myoblasts to form myotubules [50]. In addition, GAS2 can regulate p53 action and apoptosis through its inhibition of calpain [51]. Nestin is the initial intermediate filament protein expressed during myogenesis [52–54], and phosphorylation of this protein is associated with the disassembly of intermediate filaments during mitosis [53]. Nestin is also associated with increased survival of rat vascular smooth muscle cells under stress [55]. Our results are in agreement as nestin is expressed predominantly at the 18de stage in turkey skeletal muscle when myoblasts are proliferating.

Another intriguing gene was DAP, which was expressed most highly at the 18de stage and decreased throughout development for both RBC2 and F birds. The human homologue of the gene was first identified in surviving HeLa cells treated with an antisense cDNA library and the cytokine interferon-γ to induce apoptosis [56]. Very little is known about this gene other than it seemed to be a candidate for positive mediators of cell death [56, 57] and that it may regulate autophagy in injured planarians [58] and in amino acid-starved HeLa cells [59]. In addition, DAP may play a role as a regulator of turkey skeletal muscle proliferation and differentiation as shown by our group's recent work with satellite cells in vitro (Vellemen et al., manuscript in preparation).

Ca²⁺ Signaling/Muscle Function Genes

The Ca²⁺ ion acts as a second messenger in cell signaling and in skeletal muscle contraction. Aberrant Ca²⁺ signaling postmortem likely plays an important role in the occurrence of meat quality defects that can lead to significant economic losses in the pork and poultry industries. During the postmortem conversion of muscle to meat, increased Ca²⁺ release from the sarcoplasmic reticulum while the muscle temperature is still high is thought to be responsible for skeletal muscle hypermetabolism that accelerates pH decline, leading to protein denaturation and loss of protein functionality [60, 61]. A point mutation in the "halothane gene" or ryanodine receptor 1 (RYR1), the skeletal muscle sarcoplasmic reticulum Ca²⁺ release channel, has been identified in the pig [62] and has been associated with porcine stress syndrome (PSS), a sometimes fatal condition of malignant hyperthermia that is triggered by stress and is a major contributor to the development of PSE pork [60]. Genetic selection by pig breeders away from this mutation has not fully resolved the PSE problem. However, RYR1 became a logical candidate for examination in the turkey, a species that has been similarly selected for rapid lean muscle growth and also exhibits meat quality defects such as PSE. Further investigation of this protein in the turkey revealed a higher affinity for ryanodine in heavier commercial turkeys than those from the RBC2 line [63], suggesting a higher open-state probability of the Ca²⁺ release channel. Multiple alternatively spliced products were later discovered in the avian skeletal muscle ryanodine receptor isoform, αRYR, in RBC2 and commercial turkeys [64], and these alternative splicing sites are clustered in regions of the gene associated with increased frequency of mutations [65]. Recently, alternative splicing and changes in αRYR expression have also been implicated in the occurrence of PSE meat in broiler chickens [66].

Several genes involved in Ca²⁺ metabolism were identified as differentially expressed during development in this study. Calreticulin functions as a Ca²⁺ -binding chaperone and aids in maintaining Ca²⁺ homeostasis in the ER lumen of cells [67]. This protein also appears to be mainly expressed early in murine cardiac development and is essential for proper cardiac development [68]. Calmodulin is another Ca²⁺ -binding protein that sensitizes RYR1 to activation at nanomolar concentrations of Ca²⁺ and inhibits RYR1 at micromolar concentrations [69]. However, more recent work suggests that calmodulin is not an essential regulator of RYR1 [70]. Nevertheless, calmodulin plays a central role in calcium signal transduction and the developmentally related changes in gene expression suggest a key role for this protein.

Myostatin is a highly-conserved transforming growth factor-β (TGF-β) family member that strongly inhibits both hyperplasia and hypertrophy. Mutations or knockout of the gene result in an extreme increase in skeletal muscle mass in cattle and in mice [71]. In the current study, a dramatic decrease of myostatin expression was observed at 1d post-hatch, which agrees with our understanding of myostatin function as myofibers are rapidly undergoing hypertrophy at this stage. In addition, line differences at 18de and 16wk were observed, with RBC2 birds expressing myostatin significantly higher than F birds at 18de. This result could imply that F birds are able to undergo a higher degree of hyperplasia and actually produce more mature muscle fibers, leading to greater muscle mass. At 16wk, F-line birds expressed significantly more myostatin, suggesting that the faster growth rate of the heavier birds plateaus sooner than the slower-growing RBC2-line birds. Activin receptor type IIB (ActIIBR) functions as a serine/threonine kinase receptor for myostatin as well as other TGF- β family members, and binding of these ligands to ActIIBR activates the Smad signal transduction pathway to regulate gene expression [72]. Inhibition of myostatin activity by construction of a dominant-negative ActRIIB lacking a kinase domain resulted in increased muscle hyperplasia and hypertrophy in mice, as well as competition with the TGF-β family member follistatin [71]. Treatment of chicken fetal myoblasts with myostatin altered expression of genes involved in myogenic differentiation, cell architecture, energy metabolism, signal transduction and apoptosis [73], suggesting that increased myostatin expression in the current study could lead to similar changes.

Troponins are key proteins that are regulated by changes in intracellular Ca²⁺ concentrations and are responsible for striated muscle contraction. The three subunits of troponin interact with each other to bind Ca²⁺ (TNNC), inhibit myosin ATPase activity (TNNI), and to bind tropomyosin (TNNT) [74]. Transcript abundance of the fast/skeletal muscle isoform of TNNT, TNNT3, increased with development in the current study, peaking at 16wk. Interestingly, the cardiac isoform, TNNT2, was chiefly expressed in the embryonic skeletal muscle, and its expression decreased during skeletal muscle development, which is consistent with previous findings in chicken [75, 76]. Expression of another slow-twitch muscle/cardiac isoform of troponin, TNNI1, followed this pattern as well. The expression of all of the troponin isoforms is also in agreement with what is known about myofibrillar protein expression in rat skeletal muscle, as slow/cardiac forms are expressed during the development of myofibrils with a switch to fast/skeletal forms during postnatal growth and muscle regeneration [77]. These genes involved in muscle growth and Ca²⁺ signaling may be valuable candidates for changes observed in turkey skeletal muscle development.

Miscellaneous/Genes of Unknown Function

One of the goals of the current study was to identify novel genes with unknown or uncertain functions that may play important roles in skeletal muscle development. These genes may play important roles in myogenesis and development even though their annotations and functions have yet to be defined or may not have clear roles in skeletal muscle. In the current study, one gene without known annotation but identified as very highly differentially expressed by microarray analysis was chosen for further confirmation by qPCR along with others that did not easily fit into a functional category. These genes all showed line differences at 1d, indicating a possible important role during muscle hypertrophy that may help explain phenotypic differences between the slow-growing randombred line and the line selected for increased 16-wk body weight.

Several genes identified in this study were first discovered and their activities characterized in tissue types other than skeletal muscle. For example, spondin 2, (also known as mindin), which fell into the ECM group in our study, is expressed in spleen, lymph nodes, and dendritic cells and may be crucial to immune function [78]. In the current study, spondin 2 was expressed at highest levels at 1d post-hatch, over 11 times higher than its 18de expression in both turkey lines. The MGP gene was first cloned in chickens and found in highest levels in bones and relatively low levels in skeletal muscle [79]. Previous studies indicate its role in the inhibition of vascular calcification, as MGP knockout mice die within 2 months due to blood vessel rupture [80]. It may also be important in intracellular calcium homeostasis [81]. While MGP was grouped into the "Cell Death" function by the IPA software, its association with apoptosis is not clear [81]. It was another gene that appeared to be "turned on" at 1d post-hatch, with 80.28 and 55.38 times the 18de expression in RBC2 and F lines, respectively. Thus, this study has begun to uncover novel roles of these genes, which have been characterized in other tissues, in muscle growth and development.

Other interesting findings include the variation in the timing and extent of expression during development. A few stand-out genes were highly expressed at only one stage of development, clearly turned on or off in what appears to be carefully choreographed regulation. During the embryonic stage, when muscle cell hyperplasia is occurring, versican, glypican-1, betaparvin, nestin, TNNI1, and TNNT2 were expressed at least 5 times, and in the case of nestin greater than 100 times, higher than later developmental stages. Spondin 2 and MGP were both highly expressed at 1d posthatch with very little expression at 18de or 16wk. The NT5C3 and TNNT3 genes peaked at 16wk of development.