Abstract

Background:Chordoma are rare slow-growing tumors of the axial skeleton, which are thought to arise from remnants of the notochord. Little is known about the underlying mechanisms that drive this tumor. However, the assessment of gene expression levels by quantitative real-time polymerase chain reaction (qRT-PCR) is hampered due to a lack of validated reference genes. Using an unstable reference gene in qRT-PCR may lead to irreproducible results.

Methods:The expression of 12 candidate reference genes (ACTB, B2M, T, EF1a, GAPDH, HPRT, KRT8, KRT19, PGK1, RS27a, TBP, and YWHAZ) was analyzed by qRT-PCR in flash frozen chordoma samples from 18 patients. GeNorm and NormFinder algorithms were used to rank the stability of the genes.

Results:From most to least stably expressed, the top six genes found by geNorm were PGK1, YWHAZ, ACTB, HPRT, EF1A, and TBP. When analyzed by NormFinder, the top six genes were ACTB, YWHAZ, PGK1, B2M, TBP, and HPRT. GAPDH alone, which is often used as a reference gene in chordoma gene expression studies, is not stable enough for reliable results.

Conclusion:In gene expression studies of human chordomas, PGK1, ACTB, and YWHAZ are more stably expressed, and therefore, are preferred reference genes over the most often used reference gene so far, GAPDH.

Keywords: Chordoma, GeNorn, NornFinder, qPCR, reference genes

INTRODUCTION

Chordoma are rare tumors of the axial skeleton, which are thought to arise from remnants of the embryonic notochord.[ 28 30 ] Although this is a rare tumor, it usually occurs in relatively young patients, and the general prognosis is poor.[ 14 18 ] The current standard treatment is maximal surgical resection, with additional proton or carbon-ion radiotherapy.[ 1 21 31 ] Little is known about the molecular mechanisms that drive this tumor. Understanding these driving mechanisms may lead to targeted therapy which may improve the prognosis.

Quantitative real-time polymerase chain reaction (qRT-PCR) is one of the most widely used methods to study gene expression. The greatest advantages of this technique are the sensitivity and accurate quantification.[ 11 ] Because this technique is ultra-sensitive, even a small error, such as a pipetting error, variations in RNA extraction, or reverse transcription efficiency, can lead to a misinterpretation of the results.[ 3 ] Sample preparation for reverse transcriptase PCR is a multi-step process, which results in a variable RNA quality. There are several ways to evaluate the RNA quality and quantity, however, such methods provide only an indication of the total RNA quantity.[ 27 ] To correct for the errors that are accumulated by this multistep process, an internal control is usually used. This internal control (i.e., reference gene) is typically a gene that is required for the maintenance of basic cellular functions and is expressed in all cells.[ 6 ] In qRT-PCR, one or more reference genes are used as a reference point for the expression levels of other genes. By comparing the expression of a gene of interest with a reference gene in every sample, an accurate quantification of the gene of interest can be made. Ideally, the reference gene should not be regulated or influenced by the experiment. By now, multiple studies have shown that the use of several reference genes is more reliable than a single reference gene.[ 16 22 ] Thus, normalization against a single reference gene should only be applied when there is clear evidence that confirms the invariant expression under the experimental conditions.[ 4 ] Therefore, a validation of reference genes is advisable when analyzing a new tissue type or experimental condition.[ 5 ] For gene expression experiments, the choice of suboptimal reference genes may lead to irreproducible results.[ 23 ]

In chordoma, different reference genes have been used for normalization.[ 9 12 19 28 ] However, to our knowledge, none of these studies validated the stability of the applied reference genes in chordoma. Normalizing the gene expression of a gene of interest using a reference gene that is not stably expressed may lead to irreproducible results. Here, we evaluate 12 candidate reference genes in skull base chordoma samples obtained from 18 patients.

MATERIALS AND METHODS

Samples

Eighteen skull base chordoma tissue samples from 18 patients were collected in the course of resective therapy in the operation room. Upon collection, samples were flash frozen in liquid nitrogen immediately. A pathological diagnosis was made based on the World Health Organization classification of tumors. All tumors were diagnosed as a classical chordoma.

RNA isolation and cDNA synthesis

RNA was extracted from 30 mg of tissue with TRIzol^® reagent (Invitrogen, Lot no. 10531501). After phase-separation with chloroform, the aqueous phase was mixed with 1 volume of isopropanol. The RNA was pelleted and washed twice with 75%, once with 100% ethanol, and dissolved in 30 μl DEPC treated water. RNA purity and yield were estimated by optical density using a NanoDrop ND-1000 spectrophotometer (Isogen Life Science). Trace amounts of DNA were removed by DNase I with a commercially available kit (#EN0521, Thermo Scientific) using 1 μg RNA and 1 U DNase I per reaction in a reaction volume of 11 μl, according to the manufacturer's protocol. DNase-treated RNA (1 μg per reaction) was then used as a template for cDNA synthesis using a RevertAid first strand cDNA synthesis kit (Thermo Scientific, #K1622, The Netherlands) containing oligo (dT)₁₈ primers, according to the manufacturers’ protocol. A minus-reverse transcriptase (−RT) control was made for one sample to exclude genomic DNA contamination. All samples were diluted 1:5 and 1:50 in RNase free water, and then stored at −80°C until further analysis.

Primers

Twelve primers were designed to span an exon–intron boundary to exclude genomic DNA amplification. Reference genes were selected from different functional classes to avoid coregulation. The forward and reverse primer sequences are listed in Table 1 . PCR efficiency for each primer pair was measured by a standard curve via serial dilutions of cDNA. Primer sequences, efficiency, and amplicon length are listed in Table 1 .

Table 1

Primers, the calculated amplicon length and the primer efficiency for the genes analyzed in this study

qRT-PCR

Reverse transcriptase quantitative PCR (RT-qPCR) was performed in a 96-well plate using a lightcycler 480 Real-Time PCR system (Roche applied science, Mannheim, Germany). The reaction conditions were set as follows: pre-incubation (5 min at 95°C), followed by 40 cycles of 15 s at 95°C and 60 s at 60°C. Finally, a melting curve was analyzed with a ramp rate of 0.29°C/s. The reaction volume was 20 μl, with 10 μl SYBR green mastermix (Roche Diagnostics GmbH, Mannheim, Germany, #13953820), 200 nM forward/reverse primers, and 2 μl of either 1:5 or 1:50 diluted cDNA, depending on gene expression. Per reference gene, all samples were measured in the same run in duplicate. Product specificity was confirmed by analyzing by the melting peak. Two negative controls were analyzed, one containing only the reaction reagents but no template and the –RT control.

Analysis

RNA copy numbers were quantified using the comparative delta-delta-Ct method. The expression stability was tested with GeNorm[ 26 ] and NormFinder software.[ 2 ] GeNorm software ranks all tested genes with an expression stability measure (M-value). This measure is established by a pairwise comparison of the variation of one reference gene to that of the combination of all reference genes. This is based on the idea that two ideal reference genes should have an identical expression ratio in all samples, regardless of the experimental condition or cell type.[ 26 ] The lower the M-value, the more stable the gene. The geNorm algorithm also determines the optimal number of reference genes (V-value) using pairwise variation analysis. GeNorm analysis was performed with qbase plus software, version 3.0 (Licence id: 33732).

Similarly, in NormFinder, genes with the lowest stability value have the most stable expression. However, unlike the geNorm algorithm, this algorithm does not calculate an optimal number of reference genes. The Microsoft Excel-based applet NormFinder (version 0.953) was used as described in a previous study.[ 2 ]

RESULTS

Twelve potential reference genes (ACTB, B2M, T, EF1a, GAPDH, HPRT, KRT8, KRT19, PGK1, RS27a, TBP, and YWHAZ) were tested for expression stability on 18 flash frozen chordoma samples in duplicate. mRNA copy numbers were calculated using the delta–delta Ct method and the standard curve.[ 17 ] The average Ct values ranged from 22.6 in EF1a to 29.6 in TBP [ Figure 1 ]. RS27a expression was detected with the smallest standard deviation (2.3) and the most narrow range (21.7–29.1). Keratin 19 had the highest standard deviation (3.6) and the widest range (20.2–31.9).

Figure 1

Boxplot of the tested reference genes. The horizontal line marks the median, the box indicates the 25/75 percentiles and the whisker caps indicate the 10/90 percentiles. The y-axis shows the Ct values of the different reference genes

Under these experimental conditions, geNorm analysis showed that PGK1, YWHAZ, and ACTB were the most stably expressed, whereas Brachyury was the least stable gene [ Figure 2 ]. As expected for reference genes, the average M value of the 12 candidate genes was below 1.0. To select the genes that are most suitable for reference purposes, we made a pairwise variation analysis [ Figure 3 ]. By definition, geNorm applies a cut-off V-Value of 0.15. In this study, the V5/6 value is the first value below 0.15, which means there is no added value of using six reference genes over five. In order of most stable to least stable, geNorm analysis revealed that the five most stable genes are PGK1, YWHAZ, ACTB, HPRT, and EF1a.

Figure 2

GeNorm gene stability analysis. On the horizontal axis, the analyzed genes are distributed from least (left) to the most stable expressed gene (right). The stability is quantified with a GeNorm M value, as indicated on the y-axis

Figure 3

GeNorm determination of optimal number of reference genes, indicated by a GeNorm V value. This algorithm analyzes the stability when using multiple reference genes, and uses a cut-off value of 0.15 (indicated by a red line). When comparing the stability of 6 reference genes over 5 (V5/6), the value is 0.15, which means that there is no additional value of the usage of 6 reference genes. The most optimal number of reference genes is therefore 5

NormFinder analysis showed ACTB to be the most stably expressed reference gene, followed by YWHAZ, PGK1, B2M, and TBP [ Figure 4 ]. However, this program does not permit the calculation of the optimal number of reference genes necessary to obtain a reliable normalizing factor. Yet, when geNorm and NormFinder data were compared, ACTB, YWHAZ, PGK1, TBP, and HPRT were found in the top 6 most stably expressed genes of both analyses [ Table 2 ]. In addition, the 3 genes that were least stable (Brachyury, Ker8, and Ker19) were the same in both algorithms.

Figure 4

NormFinder gene stability analysis. On the horizontal axis the analyzed genes are distributed from the least (left) to the most stable gene (right). The y-axis marks the stability value

Table 2

Ranking from the most to least stable expressed top 6 reference genes as analyzed by NormFinder and geNorm

DISCUSSION

In qRT-PCR, the analysis of constitutively expressed genes is often used to calibrate the expression level of the genes of interest. As the expression of genes can vary according to tissue type and experimental condition, the choice of reference genes may warrant a validation experiment. To the best of our knowledge, there is currently no data on the required set of reference genes to study gene expression in chordoma. In this study, we evaluated 12 candidate reference genes in flash frozen chordoma samples by using the frequently used algorithms geNorm and NormFinder.

Using the cut-off V-value 0.15, geNorm analysis indicated that the optimal set of reference genes consisted of PGK1, YWHAZ, ACTB, HPRT, and EF1a. In combination with the NormFinder data, the ultimate set of reference genes is ACTB, YWHAZ, PGK1, TBP, and HPRT. These results confirm the idea that chordoma are a very heterogeneous group of tumors since not just one or two but five reference genes are required.

Glyceralde-3-phosphate-dehydrogenase (GAPDH) has been frequently used as reference gene, including in many chordoma studies.[ 9 10 12 19 ] Here, we show that GAPDH as a single reference gene is not reliable enough for qPCR analysis. Indeed, more recent studies suggest that GAPDH is not stable enough to be considered as a reference gene in many different tissues.[ 7 15 ] For instance, it is reported that GAPDH is upregulated in hypoxia,[ 8 ] and that GAPDH expression in breast cancer can be used as a prognostic marker.[ 24 ] Thus, normalizing qRT-PCR results to this reference gene may cause bias. Gene expression studies in chordoma that have used GAPDH for subtle gene expression changes, may therefore be difficult to reproduce.

In this study, we only used chordoma with brachyury expression. By immunohistochemistry, 81.1 to 100% of all chordoma are shown to be positive for brachyury.[ 13 20 28 ] It is debated whether brachyury-negative chordoma are a biologically distinct atypical subgroup.[ 13 25 28 ] In the current study, we excluded samples that did not show brachyury immunoreactivity. Even though Brachyury expression was present in all our analyzed chordoma samples, it is not stably expressed, and is thus not suitable as a reference gene. This observation is in line with a recent study in which brachyury protein levels examined by immunohistochemistry were highly variable.[ 32 ]

In this study, only tissue from skull based chordomas was used. Therefore, we cannot rule out possible differences in gene expression between cranial, sacral, or mobile spine chordoma. Another limitation of our study is that gene expression studies usually require control tissue. Since chordoma arises from remnants of the fetal notochord, this tissue can be considered to be the most suitable control. However, as notochord is difficult to dissect many studies have used nucleus pulposus tissue because it may still contain contain some notochordal cells.[ 29 ]

CONCLUSION

This is the first study that evaluated a panel of potential reference genes for chordoma. Using the geometric mean analysis of PGK1, ACTB, YWHAZ, TBP, and HPRT as reference genes allows a reliable interpretation of qRT-PCR data. The use of GAPDH as a single reference gene is not advised.

Financial support and sponsorship

Nil.

Conflicts of interest

There are no conflicts of interest.

References

1. Amichetti M, Cianchetti M, Amelio D, Enrici RM, Minniti G. Proton therapy in chordoma of the base of the skull: A systematic review. Neurosurg Rev. 2009. 32: 403-16

2. Andersen CL, Jensen JL, Orntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004. 64: 5245-50

3. Bustin SA. Absolute quantification of mRNA using real-time reverse transcription polymerase chain reaction assays. J Mol Endocrinol. 2000. 25: 169-93

4. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M. The MIQE guidelines: Minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009. 55: 611-22

5. Derks NM, Muller M, Gaszner B, Tilburg-Ouwens DT, Roubos EW, Kozicz LT. Housekeeping genes revisited: Different expressions depending on gender, brain area and stressor. Neuroscience. 2008. 156: 305-9

6. Eisenberg E, Levanon EY. Human housekeeping genes are compact. Trends Genet. 2003. 19: 362-5

7. Glare EM, Divjak M, Bailey MJ, Walters EH. beta-Actin and GAPDH housekeeping gene expression in asthmatic airways is variable and not suitable for normalising mRNA levels. Thorax. 2002. 57: 765-70

8. Higashimura Y, Nakajima Y, Yamaji R, Harada N, Shibasaki F, Nakano Y. Up-regulation of glyceraldehyde-3-phosphate dehydrogenase gene expression by HIF-1 activity depending on Sp1 in hypoxic breast cancer cells. Arch Biochem Biophys. 2011. 509: 1-8

9. Hsu W, Mohyeldin A, Shah SR, ap Rhys CM, Johnson LF, Sedora-Roman NI. Generation of chordoma cell line JHC7 and the identification of Brachyury as a novel molecular target. J Neurosurg. 2011. 115: 760-9

10. Huang B, Cohen JR, Fernando RI, Hamilton DH, Litzinger MT, Hodge JW. The embryonic transcription factor Brachyury blocks cell cycle progression and mediates tumor resistance to conventional antitumor therapies. Cell Death Dis. 2013. 4: e682-

11. Huggett J, Dheda K, Bustin S, Zumla A. Real-time RT-PCR normalisation; strategies and considerations. Genes Immun. 2005. 6: 279-284

12. Ji Z, Long H, Hu Y, Qiu X, Chen X, Li Z. Expression of MDR1, HIF-1alpha and MRP1 in sacral chordoma and chordoma cell line CM-319. J Exp Clin Cancer Res. 2010. 29: 158-

13. Kitamura Y, Sasaki H, Kimura T, Miwa T, Takahashi S, Kawase T. Molecular and clinical risk factors for recurrence of skull base chordomas: Gain on chromosome 2p, expression of brachyury, and lack of irradiation negatively correlate with patient prognosis. J Neuropathol Exp Neurol. 2013. 72: 816-23

14. Lee J, Bhatia NN, Hoang BH, Ziogas A, Zell JA. Analysis of prognostic factors for patients with chordoma with use of the California Cancer Registry. J Bone Joint Surg Am. 2012. 94: 356-63

15. Li R, Shen Y. An old method facing a new challenge: Re-visiting housekeeping proteins as internal reference control for neuroscience research. Life Sci. 2013. 92: 747-51

16. Lisowski P, Pierzchala M, Goscik J, Pareek CS, Zwierzchowski L. Evaluation of reference genes for studies of gene expression in the bovine liver, kidney, pituitary, and thyroid. J Appl Genet. 2008. 49: 367-72

17. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001. 25: 402-8

18. McMaster ML, Goldstein AM, Bromley CM, Ishibe N, Parry DM. Chordoma: Incidence and survival patterns in the United States, 1973-1995. Cancer Causes Control. 2001. 12: 1-11

19. Nelson AC, Pillay N, Henderson S, Presneau N, Tirabosco R, Halai D. An integrated functional genomics approach identifies the regulatory network directed by brachyury (T) in chordoma. J Pathol. 2012. 228: 274-85

20. Oakley GJ, Fuhrer K, Seethala RR. Brachyury, SOX-9, and podoplanin, new markers in the skull base chordoma vs chondrosarcoma differential: A tissue microarray-based comparative analysis. Mod Pathol. 2008. 21: 1461-9

21. Park L, Delaney TF, Liebsch NJ, Hornicek FJ, Goldberg S, Mankin H. Sacral chordomas: Impact of high-dose proton/photon-beam radiation therapy combined with or without surgery for primary versus recurrent tumor. Int J Radiat Oncol Biol Phys. 2006. 65: 1514-21

22. Radonic A, Thulke S, Mackay IM, Landt O, Siegert W, Nitsche A. Guideline to reference gene selection for quantitative real-time PCR. Biochem Biophys Res Commun. 2004. 313: 856-62

23. Ragni E, Vigano M, Rebulla P, Giordano R, Lazzari L. What is beyond a qRT-PCR study on mesenchymal stem cell differentiation properties: How to choose the most reliable housekeeping genes. J Cell Mol Med. 2013. 17: 168-80

24. Revillion F, Pawlowski V, Hornez L, Peyrat JP. Glyceraldehyde-3-phosphate dehydrogenase gene expression in human breast cancer. Eur J Cancer. 2000. 36: 1038-42

25. Shen J, Li CD, Yang HL, Lu J, Zou TM, Wang DL. Classic chordoma coexisting with benign notochordal cell rest demonstrating different immunohistological expression patterns of brachyury and galectin-3. J Clin Neurosci. 2011. 18: 96-9

26. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002. 3: RESEARCH0034-

27. Vermeulen J, De Preter K, Lefever S, Nuytens J, De Vloed F, Derveaux S. Measurable impact of RNA quality on gene expression results from quantitative PCR. Nucleic Acids Res. 2011. 39: e63-

28. Vujovic S, Henderson S, Presneau N, Odell E, Jacques TS, Tirabosco R. Brachyury, a crucial regulator of notochordal development, is a novel biomarker for chordomas. J Pathol. 2006. 209: 157-65

29. Weiler C, Nerlich AG, Schaaf R, Bachmeier BE, Wuertz K, Boos N. Immunohistochemical identification of notochordal markers in cells in the aging human lumbar intervertebral disc. Eur Spine J. 2010. 19: 1761-70

30. Yamaguchi T, Iwata J, Sugihara S, McCarthy EF, Karita M, Murakami H. Distinguishing benign notochordal cell tumors from vertebral chordoma. Skelet Radiol. 2008. 37: 291-9

31. Yasuda M, Bresson D, Chibbaro S, Cornelius JF, Polivka M, Feuvret L. Chordomas of the skull base and cervical spine: Clinical outcomes associated with a multimodal surgical resection combined with proton-beam radiation in 40 patients. Neurosurg Rev. 2012. 35: 171-82

32. Zhang L, Guo S, Schwab JH, Nielsen GP, Choy E, Ye S. Tissue microarray immunohistochemical detection of brachyury is not a prognostic indicator in chordoma. PLoS One. 2013. 8: e75851-