Bmc Molecular Biology Selection of Reference Genes for Expression Studies with Fish Myogenic Cell Cultures

Background: Relatively few studies have used cell culture systems to investigate gene expression and the regulation of myogenesis in fish. To produce robust data from quantitative real-time PCR mRNA levels need to be normalised using internal reference genes which have stable expression across all experimental samples. We have investigated the expression of eight candidate genes to identify suitable reference genes for use in primary myogenic cell cultures from Atlantic salmon (Salmo salar L.). The software analysis packages geNorm, Normfinder and Best keeper were used to rank genes according to their stability across 42 samples during the course of myogenic differentiation.


Background
Skeletal muscle myogenesis involves numerous steps including the proliferation, migration and fusion of myoblasts to form myotubes; the onset of myofibrillargenesis, and the maturation and hypertrophy of muscle fibres [1,2]. Myogenesis in teleost fish has several unique features compared to mammals, including the production of myotubes throughout much of adult life [3]. The in vitro culture of fish myogenic cells is an attractive system for studying the formation and differentiation of myotubes and examining the effects of various regulatory molecules on gene expression under precisely controlled conditions [4,5]. Furthermore, since traditional gene "knockouts" are unavailable in fish, cell culture provides a viable alternative for functional assays.
A pre-requisite for the quantitative measurement of gene expression is the identification of suitable reference genes to normalise the data [6,7]. Reference genes are required to normalise for differences in RNA input and mRNA/ rRNA ratios between samples [8]. Also, differences in reverse transcription efficiencies between samples can occur due to the presence of inhibitors carried over from the RNA purification [8], and the presence of PCR inhibi-tors can affect the number of cycles required to reach the quantification cycle value [9]. As gene expression patterns change in response to many stimuli, stable expression of reference genes needs to be confirmed for each experimental system. For example, genes identified as being stable in whole muscle samples, may not be suitable as reference genes in myogenic cell culture due to the vast changes in cell metabolism and structure that occur during the transition from myoblast to myotube. Previous myogenic cell culture experiments using the C2C12 cell line have relied on Actb [10] and Gapd [11] as internal reference genes, however, the validity of these genes is questionable as Gapd [12,13] and Actb [14] expression has been shown to vary considerably.
A number of computer based analysis packages have been developed which analyse gene expression patterns and allow for the identification of stable reference genes. Vandesomple et al [15] designed the widely used geNorm package which uses a pairwise analysis of gene expression to identify stable reference genes. Likewise, Bestkeeper [16] performs a pairwise comparison, whereas Normfinder [17] uses a mathematical model to estimate overall expression variation of candidate reference genes, but also the variation between sample groups. Vandesomple et al [15] demonstrated that use of a single reference gene can lead to aberrant gene expression values, and now it is widely accepted that using several reference genes for normalisation is preferable.
Currently there is no information available on reference gene stability in fish myogenic cell cultures. In this paper we examine the stability of eight potential reference genes during the transition from single nucleated myoblasts to multinucleated myotubes in myogenic cell cultures derived from Atlantic salmon, one of the most commercially important aquaculture species.

Results
Cell cultures were visualised using confocal microscopy and the phenotype of cells determined at 2 d, 5 d, 8 d, 11 d and 14 d (Figure 1). The myogenic nature of the cell culture was confirmed by the presence of the myogenic marker desmin (Figure 1 a, e, i, m, q) and the presence of multi-nucleated myotubes visualised by Alexa Fluor 568phalloidin stained actin filaments (Figure 1 b, f, j, n, r) and nuclei stained with sytox green (Figure 1  Each of the candidate reference genes tested gave amplification from cDNA derived at each time point of the cell culture, while the no template control (NTC) and minus reverse transcription controls (-RT) gave no signal. The specificity of each primer was by confirmed by the presence of a single band on agarose gel electrophoresis and the presence of a single peak in the dissociation curve analysis which exactly matched the dissociation curve of a plasmid standard of known sequence. Amplification of the correct product was confirmed in each case through the sequence analysis of cloned PCR products.

Reference gene stability
The cell culture undergoes many structural and metabolic changes during the transition from mononucleic cells to multi-nucleated myotubes. We therefore chose to analyse the expression of genes from early time points as well as time points after the culture has produced multi-nucleated myotubes (17 d and 20 d). The analysis of reference gene stability can thus be performed in three phases, the first in developing myotubes (2 d -11 d), the second in established myotubes (11 d -20 d) and the third covering all time points. Based on the raw expression data ( Figure  2), Sdha and Pgk show higher levels of variance than the remaining genes and appear the least stable. Figure 3 shows the raw expression values obtained for each gene at each of the time points sampled. For several of the genes, there is higher intergroup variation indicating that these genes are differentially regulated during the progression of the cell culture. For example, 18SrRNA, Pgk, and Actb all have higher Cq values at 2 d and 5 d when the culture is predominantly mononucleic cells than they do at later time points when myotubes have formed. Sdha has high inter-assay and intra-assay variation and is clearly unsuitable as a reference gene in Atlantic salmon myogenic culture.
As inspection of raw Cq values alone is insufficient for determining gene expression stability, the data obtained were further analysed using three software packages Bestkeeper, geNorm and Normfinder. Each package uses a different algorithm to determine the most stable reference gene, and as no single method has been accepted as the most appropriate for identifying stable gene expression, all three packages were used for analysis.

geNorm Analysis
Data analysis using geNorm was performed two ways. The first method used the absolute values derived from a plasmid standard curve as input, the second used the delta Cq method, with the PCR efficiencies based on a dilution series of pooled cDNA samples. The results from the three geNorm analyses covering all time points (absolute value method), developing myotubes and established myotubes are shown in figure 4. When all samples were analysed, the genes were ranked in an identical order using both analysis methods from most to least stable: Growth and differentiation of myogenic cells extracted from Salmo salar fast myotomal muscle

Hprt1>RNApolII>Ppia>Ef1α>18SrRNA>Actb>Pgk>Sdha.
Analysis of days 2-11 using both analysis methods revealed the same order of stability as when all days were analysed except for Ef1 and Ppia swapping order. When days 11-20 are analysed using the absolute method, the order changes from most to least stable: Pgk>Actb>Hprt1>Ppia>Ef1>18SrRNA>RNApolII>Sdha. When the delta Cq method was used, the order changed to:Pgk>Hprt1> Actb>Ppia>Ef1>18SrRNA>RNApolII>Sdha. The change in order from days 2-11 to days 11-20 likely reflects the changes in metabolism and structure that occur during differentiation and growth process as myotubes form. In all three analyses, the M values obtained for Hprt1, Ppia and Ef1α were quite similar, ranging from 0.23-0.36. Based on the similar M values, it would appear that any combination of Hprt1, Ppia and Ef1α would be suitable for normalisation.

Bestkeeper Analysis
Using the initial statistics produced by Bestkeeper ( Figure  2), the genes were ranked in the following order from most to least stable: Hprt1>RNApolII >Ppia>Ef1α>Sdha> Actb>18SrRNA>Pgk when all time points were examined. For 2d-11d, the genes were ranked: Hprt1>RNApolII>Ef1α>Ppia>Actb>18SrRNA>Sdha>Pgk and for 11d-20d: Pgk>Hprt1>Ppia>RNApolII> Ef1α>Actb>18SrRNA>Sdha. All candidate reference genes examined were positively correlated with each other (Table 1), with the highest correlations found between Actb/Pgk (r = 0.914) and Ppia/Actb (r = 0.874). Correlations between the remaining genes ranged from 0.198 for RNApolII/Ef1α to 0.814 for 18SrRNA/Actb ( Table 1). The low level of correlation between many of the genes is due to the small inter and intra-group variation observed for the majority of the genes. From the initial statistics, the four least stable genes were removed from further analysis. The algorithm used in Bestkeeper then calculates the correlation of each gene with the Bestkeeper Index which is calculated as the geometric mean of the candidate reference genes. The candidate reference genes were ranked in order based on their correlation with the Bestkeeper index based on the four most stable genes. For all time points the genes were ranked in order most to least stable: Ppia>Ef1α>Hprt1>RNApolII whereas in developing myotubes the order was: Ppia>Ef1>RNApolII>Hprt1 and in established myotubes the ranking was Ef1α>Ppia >Pgk>Hprt1.

Normalisation
Based on the results from the three analysis methods, four genes, EF1α, Ppia, Hprt1 and RNApolII are consistently stable. In order to assess the stability of the normalisation factors obtained, we first compared the normalised expression of Des to various combinations of the geometric average of two genes ( Figure 5). Six normalisation factors were derived by calculating the geometric averages of the following gene combinations: A: HPRT1, PPIA; B: RNApolII, HPRT1; C: EF1α, HPRT1; D: EF1α, Ppia; E: RNApolII, EF1α; F: RNApolII, Ppia. We found that there were significant differences (ANOVA P < 0.05) in Des gene expression at day 11 (B v D), at day 17 (A v B, D, and F; B v C, D and E; C v E and F; D v E and F) and at day 20 (B v C, D and E; D v F). We therefore examined normalisation using geometric average of three genes ( Figure 6). Four normalisation factors were derived by calculating the geometric averages of the following gene combinations: A:   that using the geometric average of any three of these genes is suitable for normalisation.

Discussion
In this study, we examined the expression of eight candidate reference genes for normalisation of quantitative real-time PCR data from a primary culture of Atlantic salmon myogenic cells. The identification of genes with stable expression in all samples of an experiment is crucial as it is necessary to normalise for variability between samples introduced during the production of the cDNA [6,7]. As a universal reference gene with stable expression in all experimental systems is not available, suitable reference genes for each experiment need to be determined.
Myogenic cell culture is characterised by distinct phases where cells first proliferate, and then fuse to form multinucleated myotubes [2]. We therefore identified genes that were stable early time points where the majority of the cells are mononucleic and forming small myotubes (culture days, 2-11) those stable in established myotubes (days [11][12][13][14][15][16][17][18][19][20] and those most stable for the entire culture period. When the raw Cq values obtained at each time point are compared (Figure 3), it is clear that for the majority of the genes examined, stable expression is observed once myotubes have become established in the culture, whereas higher intra and intergroup variation is observed when myotubes are developing. For example, Pgk, a glycolytic enzyme, appears to be upregulated as  myotubes start to form, and then has stable expression in established myotubes as indicated by the low inter and intra-group variation (Figure 3).
To identify stably expressed genes, analysis packages such as geNorm and Bestkeeper perform a pairwise comparison of gene expression across the various samples in an experiment. Therefore it is crucial that genes used are not coregulated or present on the same pathway, as co-regulated genes will likely have similar expression patterns and would therefore appear to be stably expressed in any biological experiment. For this reason we chose genes involved in a number of different biological processes (Table 2), such as nucleotide recycling (Hprt1), peptide isomerisation (Ppia), glycolysis (Pgk), citric acid cycle (Sdha), ribosome assembly (18SrRNA), transcription (RNApolII), translation (Ef1α) and cytoskeleton structure (Actb).
Based on the M values obtained in geNorm (Figure 4), the stability index from Normfinder ( Figure 4) and the descriptive statistics produced by Bestkeeper (Figure 2), it would appear that several of the genes used in this study are suitable for normalisation of gene expression data from salmon myogenic cell cultures. For example, Pfaffl et al [16] recommends using genes that have a standard deviation for the Cq values less than one for calculating a Bestkeeper index. In our study, all genes examined had standard deviation less than one. This is also reflected in the slight changes in the order of gene stability obtained from each of the three software packages. The least stable gene identified by all analysis methods was Sdha. Sdha has been used as a reference gene in a number of studies using different tissues [18,19], however its high inter and intragroup variation make it unsuitable for normalisation in salmon myogenic cell cultures.
Results obtained from geNorm identify Hprt1 and RNApolII as the most stable genes when all time points were examined, however, the M values obtained for Ppia and Ef1α are quite similar and thus these genes are also likely to be suitable for normalisation. The same set of genes was found to be stable in developing myotubes, but differed in established myotubes where Pgk/Actb were found to be the most stable, although the Hprt1 and Ppia also had low M values and can be considered stably expressed.
The most stable genes identified for all time points by Normfinder ranked in descending order were Ppia>18SrRNA>Hprt1>Ef1α. Both Ppia and Hprt1 have been reported to give stable expression in mouse C2C12 Normalisation of desmin mRNA expression to various com-binations of the geometric average for two genes   myotubes [20,21] and Ef1α has been reported to have stable expression in some Salmon tissues [22]. As the 18S and 28S ribosomal RNAs are highly abundant and account for the vast majority of RNA, it is unsurprising that 18SrRNA is found to be stable across the samples as equal amounts of RNA were reverse transcribed. However, Vandesomple et al [15] criticise the use of 18SrRNA as a housekeeping gene due to its high abundance making baseline subtraction difficult. Also, transcription of rRNA and mRNA occur via RNA polymerase I and II respectively which may lead to imbalances in the two mRNA fractions as reported by Solanas et al [23]. The similar stability indices obtained for Ppia, Hprt1 and Ef1α, identify all of these genes as suitable for normalisation.
Similar to the results of geNorm and Normfinder, Bestkeeper analysis revealed the most stable genes to be Ppia, Ef1α, Hprt1 and RNApolII. Interestingly, Actb, which has been used as a reference gene in numerous studies [10,24], was found to be the third least stable gene in this analysis, having high intra-group variation in developing myotubes as well as high inter-group variation when comparing developing and established myotubes. These differences in Cq values between developing and established myotubes indicate that Actb is differentially regulated during differentiation of Atlantic salmon myogenic cells, as reported in chicken and mouse myoblast culture [25,20] and is therefore unsuitable as reference gene for myogenic culture. Interestingly, RNApolII and Hprt1, which were identified as the most stable genes in geNorm ( Figure 4) had a correlation coefficient of only 0.57, which was lower then for many of the other genes ( Table 1). The selection as the most stable genes in geNorm is likely a reflection of the low intra and inter-group variation observed for both of these genes (Figure 3).
Vandesomple et al [15] recommend using the geometric average of three reference genes for accurate normalisation. To assess the suitability of the reference gene candidates, we first normalised the expression of Des to combinations of the geometric average of two reference genes from Ppia, Hprt1, Ef1α and RNApolII ( Figure 5). We found significant differences in Des expression at days 11, 17 and 20 when comparing results from different combinations of reference genes. However, when three genes were used, there were no significant differences between any of the combinations of reference genes ( Figure 6) indicating that all four genes are suitable for normalisation when the geometric average of three genes is used.

Conclusion
To the best of our knowledge, this is the first study examining gene expression stability in myogenic culture of a teleost species and thus provides a useful platform for gene expression studies using this system. The data provided in this paper may also be useful in guiding researchers performing myogenic cell culture in other teleost species. We recommend using a three gene normalisation factor using the geometric average of any combination of EF1α, Ppia, RNApolII and Hprt1.

Isolation of myogenic satellite cells
Myosatellite cells were isolated using a method similar to that described by Koumans et al [26]. Juvenile Atlantic salmon (Salmo salar L) 30 ± 6 g (mean ± s.d., N = 10) were used for each culture. As the experimental animals had not undergone gonadal development, the gender of the fish was not determined. Fast myotomal muscle was dissected under sterile conditions and placed in extraction media consisting of Dulbecco's modified eagle's media (DMEM) 9 mM NaHCO3, 20 mM HEPES (pH 7.4) with 15% (v/v) horse serum and 1 × antibiotics (100 units/ml penicillin G, 100 μg/ml streptomycin sulfate, 0.25 μg/ml amphotericin B) (Sigma, Gillingham, Dorset, UK) at a ratio of 1 gram of muscle per 5 ml extraction media. The tissue was then minced with a sterile scalpel before centrifugation at 300 g for 5 min, and two washes with DMEM without horse serum. The muscle pieces were digested with collagenase (0.2% m/v in DMEM, Type 1a, Sigma, Gillingham, Dorset, UK) for 70 minutes at room temperature in the dark, before centrifugation at 300 g for 5 minutes. The resulting pellet was washed twice with DMEM before being passed through a pipette repeatedly to separate cells. Samples were further digested with trypsin (0.1% in DMEM) for 20 minutes at room temperature. The resulting cell suspension was centrifuged (300 g, 1 min). The supernatant was poured into 20 × vol of extraction media containing serum to inhibit trypsin activity. The pellet was further digested by a second treatment with trypsin for 20 min at room temperature, before centrifugation at 300 g 1 min. The supernatant was poured into 20 × volume of extraction media. The extraction media containing the cell suspension was centrifuged 300 g, 20 min. Cell pellets were re-suspended in 30 ml of basal medium before mechanical trituration through 10 ml and 5 ml pipettes until cells are separated. The cell suspension was then passed through 100 μm and 40 μm nylon cell strainers (BD Biosciences San Jose, CA, USA) and centrifuged 20 min 300 g. The cells were resuspended in basal media, cell number determined using hymaecytometer, and then diluted to give approximately 1.5 × 10 6 cells/ml.

Cell culture
All cell culture methods were performed using Aseptic technique in a Microflow 2 Advanced biosafety cabinet (Bioquell Ltd, Andover, UK). 6 well cell culture plates (Greiner Bio-One Ltd, Stroudwater, UK) were treated with a 100 ug/ml poly-lysine solution (Sigma, Gillingham, Dorset, UK) at 4 μg/cm 2 for 5 minutes at room temperature, then aspirated before 2 washes with sterile water and allowed to air dry. 1 ml of laminin (Sigma, Gillingham, Dorset, UK) in DMEM at 20 μg/ml was applied to each well and incubated at 18°C overnight prior to plated cells. Cell culture was performed using complete medium (DMEM, 9 mM NaHCO3, 20 mM HEPES (pH 7.4), supplemented with 10% foetal calf serum (Sigma, Gillingham, Dorset, UK) and 1 × antibiotics (Sigma, Gillingham, Dorset, UK) which was changed daily.

Immunofluorescence of culture cells
Cells were grown on glass coverslips treated with poly-Llysine and laminin as described above.

Quantitative real time PCR experiments
The following procedures were performed as to comply with the MIQE guidelines [27].  Hamburg, Germany). The concentration of each plasmid was calculated based on absorbance at 260 nm, and a dilution series produced for calculation of copy number via qPCR.

Primer design
Primers were designed using NetPrimer (Premier BioSoft, Palo Alto, CA, USA) to have Tm of 60°C, and where possible, were designed to cross an exon-exon junction to avoid amplification of contaminating genomic DNA. To determine exon-intron junction sites, genomic sequences for orthologous genes from Danio rario, Gasterosteus aculeatus, Oryzias latipes, Takifugu rubripes and Tetraodon nigroviridis were retrieved from Ensembl http:// www.ensembl.org/index.html, and compared to the Salmo salar cDNA sequences using the Spidey software tool http://www.ncbi.nlm.nih.gov/spidey/. Primers were designed across conserved exon-intron junctions and used at a final concentration of 500 nM. 18SrRNA, at 500 nM gave poor amplification efficiencies, however this was improved using a final concentration of 1.5 μM. The primers used for qPCR are listed in table 3 and have been submitted to rtprimerdb http://www.rtprimerdb.org/ [28].

Data analysis
The stability of candidate reference genes was determined using geNorm [15], Normfinder [17] and Bestkeeper [16]. Input data for geNorm and Normfinder were absolute values derived from a plasmid standard curve with the data for geNorm transformed as per author's guidelines. Input for Bestkeeper was the Cq values, and the PCR efficiencies calculated from a dilution series (1/20, 1/40, 1/80, 1/160, 1/320, 1/640) of cDNA. Normfinder Analysis of inter and intra group variation was performed on all data, days 2-11 and days 11-20. Statistical analysis was performed with Minitab (Minitab Inc).