Regular paper

STUDY OBJECTIVE
The aim of this study was to test a panel of 6 reference genes in order to identify and validate the most suitable reference genes for expression studies in paired healthy and non-small cell lung cancer tissues.


METHOD
Quantitative real-time PCR followed by the NormFinder- and geNorm-based analysis was employed. The study involved 21 non-small cell lung cancer patients.


RESULTS
The analysis of experimental data revealed HPRT1 as the most stable gene followed by RPLP0 and ESD. In contrast, GAPDH was found to be the least stable gene. HPRT1 together with ESD was revealed as the pair of genes introducing the least systematic error into data normalization. Validation by bootstrap random sampling technique and by normalizing exemplary gene expression data confirmed the results.


CONCLUSION
Although HPRT1 and ESD may by recommended for data normalization in gene expression studies on non-small cell lung cancer, the suitability of selected reference genes must be unconditionally validated prior to each study.


INTROduCTION
Despite enormous effort and significant progress in therapy, non-small cell lung cancer (NSCLC) still remains a highly lethal disease with only a poor prognosis: the 5-year survival rates remain lower than 15% (Tsuboi et al., 2007).In the USA population, the mortality rate due to lung cancer was predicted to reach 30% by 2007 (Jemal et al., 2007).To increase the survival rates of NSCLC patients, the disease must be detected as early as possible.Quantitative real-time RT-PCR (qPCR) was recently shown to be useful for early NSCLC diagnosis, prognosis prediction and gene expression analysis (Hayes et al., 2006;Ford et al., 2007;Su et al., 2005;Wang et al., 2005;Zhang et al., 2005a;Cheng et al., 2005;Plummer et al., 2005;Buttitta et al., 2005;Dworakowska et al., 2004;Sozzi et al., 2003).In gene expression studies, however, to accurately and reliably quantify the target gene expression level, qPCR data need to be normalized in order to get rid of the unspecific variability arising from the differences in cDNA quality and/or quantity, using one or more reference gene(s) amplified simultaneously.Nowadays it seems, unfortunately, that an ideal reference gene, showing the same expression level irrespective of sample treatment and experimental design, either does not exist, or has not yet been discovered (Vandesompele et al., 2002;Andersen et al., 2004;Dheda et al., 2004;Radonic et al., 2004).As a consequence, the identification of reliable and stable reference gene(s) is a prerequi-

2009
P. Gresner and others site prior to any reliable interpretation of the qPCR data.
Several studies have already identified and validated suitable reference genes for gene expression studies in various human tissue and cell types (Gabrielsson et al., 2005;Ohl et al., 2005;2006;Zhang et al., 2005b;de Brouwer et al., 2006;Silver et al., 2006;Jung et al., 2007).In this study, we built a panel of six reference genes comprising four commercially available and two putative reference genes (Table 1) and analyzed them with regard to their suitability as reference genes for expression studies in paired healthy and tumor tissue samples from patients with NSCLC.In order to define the most stable reference genes the expression changes between and within these two groups were investigated and analyzed using two free software products available online and subsequently validated by a bootstrap random sampling technique.

MATeRIAl ANd MeThOdS
Tissue samples.Primary tumor and corresponding healthy lung tissues were obtained from NSCLC patients (n = 21; 9 women, 12 men) aged 50-79 years (mean age 64.69 yrs, S.D. 7.95 yrs) undergoing curative resection surgery with adjuvant chemotherapy between May, 2006 andApril, 2007.Surgically removed tumor samples, including 9 adenocarcinomas (AC; 2 WHO grade I, 7 WHO grade II), 12 squamous cell carcinomas (SqCC; 10 WHO grade II, 2 WHO grade III) together with paired healthy lung tissue samples, were cut into small pieces (approx.30 mg each), placed immediately in RNAlater RNA stabilization reagent (Qiagen, Hilden, Germany), processed according to the manufacturer's instruction and finally stored at -20°C until processing.Pri-or written informed consent for participation in the study was obtained from all individuals enrolled.The study was performed under the guidelines of the Helsinki Declaration for human research and approved by the Local Bioethics Committee for Scientific Research (resolution no. 2/2006).
The primer sequences used for expression assays (Table 2) were designed using the Primer3 primer design software (Rozen & Skaletsky, 2000).Sequence variability among individual subjects was avoided by selecting sequence regions free of known polymorphisms.The primers were designed to amplify all splicing variants.The primer specificity was confirmed by 2% agarose gel electrophoresis as well as by the melting curve analysis.
Additionally, the human glutathione peroxidase 1 (hGPX1) gene was amplified in this study along with the candidate reference genes in order to demonstrate the effect of reference gene selection on data normalization and result interpretation.The primer sequences used for hGPX1 amplification also presented in Table 2.
RNA extraction and cdNA synthesis.Total RNA was isolated from 20-25 mg of paired healthy and tumor lung tissue specimens using the RNeasy

Reference genes in lung cancer
Plus Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instruction.Genomic DNA contamination was removed by the on-column digestion with the RNAse-free DNase set (Qiagen, Hilden, Germany).Total RNA was further quantified and analyzed with regard to protein content using an Eppendorf BioPhotometer (Eppendorf, Germany) instrument.All samples with 260 nm-to-280 nm absorbance ratio higher than 1.75 were stored at -80 o C and utilized within 1 month.Purified and frozen RNA (300 ng) was then reverse-transcribed in a 30 μl mixture using an Omniscript RT Kit (Qiagen, Hilden, Germany) and Oligo dT Primers (Promega, Madison, WI, USA) according to the manufacturer's instruction on a BioRad's PTC-200 DNA Engine thermal cycler (BioRad, Hercules, CA, USA).Finally, the reverse transcription was stopped by incubating the reaction mixture at 95 o C for 5 min.
Real-time PCR.The expression levels of six reference genes (ESD, GAPDH, HPRT1, POLR2A, RPLP0 and YAP1) were measured by quantitative real-time PCR (qPCR) on an iQ5 iCycler Multicolor Real Time PCR detection system (BioRad, Hercules, CA, USA) utilizing an IQ SYBR Green Supermix (BioRad, Hercules, CA, USA).Samples were amplified in 20-μl aliquots containing 25 ng cDNA, 1× iQ SYBR Green Supermix and 100 pmol of forward and reverse primer each.The cycling conditions comprised 3 min polymerase activation at 95 o C followed by 40 cycles of denaturation at 95 o C for 15 s, annealing at 57 o C for 45 s and extension at 72 o C for 1 min.All samples were measured in triplicate and paired tumor and healthy lung tissue specimens were always analyzed in one analytical run in order to avoid between-run variations.Each analytical run was accompanied by a no-template control for each analyzed gene as well as a sample containing cDNA pooled from all experimental samples serving for a between-run control.The between-run variation was assessed by the variation coefficient of the threshold cycle (CT) for HPRT1 in pooled sample (in triplicate for each run), and it did not exceed 1.0%.data analysis.The threshold cycle (CT) and the PCR reaction efficiency were calculated for each sample based on its amplification curve using a PCR DataMiner algorithm (Zhao & Fernald, 2005) and averaged PCR efficiency for each gene was further confirmed by a dilution-series experiment (not shown).In all subsequent analyses, the average CT value from each sample's PCR triplicates and average efficiency calculated from all samples for each gene were used.Fold changes in expression quantities between the tumor and healthy lung tissue samples were calculated according to the delta-CT method based on the formula R = (1 + E) -ΔCT , where ΔCT stands for the difference between fractional threshold cycle numbers in paired tumor and healthy lung tissue samples, respectively, and E is the averaged PCR efficiency for a given gene (Schmittgen & Zakrajsek, 2000;Dheda et al., 2004).
Prior to the assessment of the candidate reference genes expression stability, Levene's test for equality of variance and paired t-test of equality of means were used to assess the equality between CT values of tumor and normal groups of lung tissues.Spearman's rank correlation was used to assess correlations between expression levels of candidate reference genes.Finally, in order to analyze both the inter-group and intra-group expression variation of the six candidate reference genes, the data were submitted to NormFinder (Andersen et al., 2004).The combined effect of both inter-and intra-group expression variations was then expressed by means of "stability value"; the lower the "stability value" the higher the expression stability of a given gene.Additionally, the candidate reference gene expression stability was also analyzed by geNorm algorithm (Vandesompele et al., 2002).data validation.We performed a bootstrap analysis in order to assess the certainty of the candidate reference gene ranking.The ranking method was bootstrapped utilizing the technique of resampling with replacement from the original set of 14 samples (Boos, 2003).The bootstrap procedure was performed 500 times.
Additionally, the effect of the reference gene used for qPCR data normalization on resulting expression level of the gene of interest and interpretation of such a result was demonstrated utilizing the mRNA expression data obtained for human glutathione peroxidase 1 (hGPX1) in healthy lung and non-small cell lung cancer tumor tissue.Following hGPX1 amplification, the normalized relative expression level (NRQ) of hGPX1 in tumor versus paired healthy sample was calculated based on the method described by Pfaffl and coworkers (2002).For qPCR data normalization, either the worst or the best combination of reference gene(s) according to the ranking revealed by NofmFinder in our study was used.

Candidate reference gene expression levels and ranges
In general, the six candidate reference genes showed a wide range of expression levels with CT values falling between 16 and 30 cycles.The median CT values were in the range usually covered by reference genes, varying from 18.2 (RPLP0) to 26.3 (POLR2A) and the genes could be divided into two groups: a group of genes with median CT value below (RPLPO and GAPDH) or above 20 (ESD, HPRT1, POLR2A and YAP1; Fig. 1).There was a trend of higher median CT values of RPLP0, HPRT1 and GAPDH in healthy lung tissue compared to paired cancer tissue.Although the difference was not big enough to assume the inequality of the expression levels between healthy and tumor tissue for RPLP0 and HPRT1, the equality of expression level of GAPDH in healthy and tumor tissue cannot be assumed (P < 0.01).For POLR2A, YAP1 and ESD, there was a trend of lower median CT values in healthy compared to tumor lung tissue, but the difference was not big enough to assume inequality of expression levels of these genes between normal healthy and tumor lung tissue.Individual candidate reference genes also showed different expression ranges across all studied samples.ESD showed the smallest variability of CT values in both healthy (min-max range: 1.0) and tumor (min-max range: 1.9) lung tissue, while YAP1 (min-max range: 5.0) and GAPDH (min-max range: 3.6) showed the largest variance in healthy normal and tumor lung tissue, respectively.

Candidate reference genes expression stability
A NormFinder algorithm was used to rank the six lung tissue candidate reference genes according to their expression stability (Table 3).As revealed by the analysis, HPRT1, characterized by the stability value of 0.186, showed the most stable expression level and thus was depicted by the algorithm as the best choice for a single reference gene for expression studies on healthy and tumor lung tissue.RPLP0 and ESD were ranked the second and third best candidate reference genes, respectively, while one of the most commonly used known reference genes, GAPDH, was ranked last following YAP1 and POLR2A, suggesting that it cannot be considered a reliable reference gene for expression assay studies on lung tissue.The program additionally identified HPRT1 and ESD as the best combination of two reference genes, with an improved stability value of 0.177 as compared to the one of the best gene alone.
These data were supported by the comparison of median fold change of expression levels between healthy and tumor lung tissue (Fig. 2).RPLP0 and HPRT1 were characterized by the lowest median Reference genes in lung cancer fold change (1.3 and 1.7, respectively), while GAP-DH showed the second highest median fold change (> 2.4).The highest median fold change of expression level between healthy and tumor lung tissue was found in POLR2A (> 3.4).Concerning the fold change variability, expressed by means of the maximal fold change observed among the samples, maximal values were recorded for POLR2A and YAP1 (> 13.2 and 12.7, respectively).On the other hand, the smallest variability was revealed for RPLP0 and ESD (maximal fold changes < 3.6 and 3.9, respectively).
Stability of the expression level of candidate genes in paired healthy and tumor lung tissue was also investigated by geNorm.The curve presented on Fig. 3A represents the stepwise exclusion of the least stable candidate reference gene.This algorithm confirmed the results obtained by NormFinder, excluding GAPDH and POLR2A at first as the least stable reference genes while identifying RPLP0 and HPRT1 as the most stable gene pair (M = 0.488).
The geNorm analysis also allows for evaluation of the optimal number of reference genes required for a reliable and accurate normalization of expression data.This number is determined based on the average pairwise variation V n /V n+1 between two sequential normalization factors (NF n /NF n+1 ) calculated from n and n+1 reference genes.Vandesompele et al. (2002) suggest that the V n /V n+1 cut-off value of 0.15 should be considered as a limit beneath which the involvement of additional reference genes would not be required.However, based on our experimental data, the variation V n /V n+1 did not reach this cut-off; observed variability ranged from 0.272 for addition of third reference gene (n = 2) to 0.175 for addition of sixth reference gene (n = 5) (Fig. 3B).

Mutual correlation among candidate reference genes
Despite the wide range of expression levels observed for the analyzed candidate reference genes in our study, we found significant correlations between the expression levels of some of the analyzed genes (Table 4).The expression levels of HPRT1, GAPDH and RPLP0 were significantly correlated with each other (P < 0.05).The expression level of YAP1 was also significantly correlated with that of POLR2A (P < 0.0001) and ESD (P < 0.05).

Validation
The bootstrapping of the ranking method revealed that the HPRT1 was ranked first in 53% of 500 bootstrap samples, followed by RPLP0 and ESD, which were ranked first in 32% and 25% of bootstrap samples, respectively (Fig. 4).GAPDH was confirmed to be the last in candidate reference gene ranking being ranked sixth in 76% of bootstrap samples, preceded by POLR2A (ranked fifth in 49% of  In order to demonstrate the direct effect of reference gene selection on qPCR data normalization, hGX1 mRNA in tumor and normal lung tissue (median CT values: 19.2 (tumor) vs. 18.6 (normal)) was amplified and the expression level in tumor versus normal healthy lung was normalized using either HPRT1 and ESD (the best combination) or GAPDH (the worst candidate gene) used as a reference.Using HPRT1 and ESD as a reference, the analysis revealed that the hGPX1 mRNA expression, despite being lowered in tumor tissue, did not vary significantly from the one observed in healthy lung tissue (NRQ = 0.589; 95% CI: 0.126-1.650;NS).On the other hand, using GAPDH for qPCR data normalization, the expression of hGPX1 was shown to be significantly down-regulated in tumor lung tissue compared to healthy lung (NRQ = 0.258; 95% CI: 0.060-0.787;P < 0.05).

dISCuSSION
In our present paper we made an attempt to identify genes with the lowest expression variability and validate their use for qPCR data normalization in expression studies on NSCLC tissue.In order to   The empirically obtained distributions of ranks shown in percentage of the candidate reference genes as revealed by bootstrapping of the ranking method utilizing the technique of resampling with replacement (size 500).Stable genes show low ranks.Data are given as percentage of bootstrap samples in which a given rank for a given gene was recorded.

Reference genes in lung cancer
normalize such data, the use of an internal RNA template has been the most commonly used approach although the choice of a suitable reference gene for a given gene expression study might pose a fundamental problem (Bustin, 2002).Nowadays, the use of a minimum of two or three housekeeping genes as a reference for qPCR data normalization is believed to be a golden standard allowing for a reliable and meaningful interpretation of findings.Moreover, it is necessary to define the scientific problem to be investigated prior to the identification of a suitable reference gene(s), since differential gene expression studies investigating the effect of various experimental and physiological conditions on gene expression profiles should consider that not only genes of interest but also reference genes could be affected (Ohl et al., 2005).We found HPRT1, RPLP0 and ESD to be the most stably expressed genes and thus identified as the most suitable reference genes for qPCR data normalization in healthy and tumor lung tissue.On the other hand, GAPDH, one of the genes which are most commonly used as endogenous controls, was ranked last.In literature, there is an increasing number of reports documenting the variable mRNA levels of GAPDH under various physiological and pathological circumstances (Zhong & Simons, 1999;Ke et al., 2000;Suzuki et al., 2000;Revillion et al., 2000;Tatton et al., 2000;Zhu et al., 2001;Goidin et al., 2001;Bustin, 2002;Bas et al., 2004;Johansson et al., 2007), nevertheless, this gene still continues to be the reference gene of choice in many recent gene expression studies without its previous validation.However, some other studies showed that following validation, this gene may be successfully used for qPCR data normalization in some instances (Winer et al., 1999;Bustin, 2002;Wall & Edwards, 2002;Meldgaard et al., 2006).
We employed two techniques in order to validate the results we obtained.While the bootstrap randomization technique confirmed the robustness of the reference gene ranking obtained by both geNorm and NormFinder, the qPCR data normalization using the worst (GAPDH) and the best combination (HPRT1 and ESD) of reference gene(s) revealed what far-reaching impact on interpretation of experimental results an incorrect reference gene selection may have.Based on qPCR data normalized by HPRT1 and ESD, one may conclude, that there are no differences in hGPX1 expression between healthy and NSCLC lung tissue, unlike the situation when GAP-DH is used as a reference.Here, the results would indicate that the hGPX1 mRNA is downregulated in tumor compared to healthy lung tissue and would thus lead to a completely different conclusion.This is of utmost interest since such incorrect qPCR data normalization may generate misleading conclusions concerning the pathophysiology of a given disease and lead scientists in the wrong direction with respect to further research on a given topic.
Up to date, two studies addressing the problem of identification of suitable reference genes for lung cancer gene expression studies have been published, but their outcomes are inconsistent (Liu et al., 2005;Saviozzi et al., 2006).Of note might be the fact that both these studies included GAPDH and RPLP0, but they were ranked differently: while in one study they were ranked among the best three genes (Liu et al., 2005;Saviozzi et al., 2006), the other revealed them as absolutely unsuitable for use as normalizers of qPCR data (Liu et al., 2005;Saviozzi et al., 2006).In the latter study the authors have suggested that the discrepancies may have arisen due to the inappropriate statistical analysis used for the assessment of the reference gene ranking in the study by Liu and coworkers (Liu et al., 2005;Saviozzi et al., 2006).Considering the low and high rankings of GAPDH and HPRT1, respectively, found in our study, our findings seem to partially agree with the results from both these earlier studies.However, it might be possible that, besides the statistical approach employed, there are some other factors contributing to the discrepancies among these three studies, including ours.First of all, the different outcomes of the three studies may be accounted for by an interindividual variability in gene expression profiles of subjects involved in each study.Second, it is possible that SqCC and AC samples would need to be treated separately in order to establish the best reference gene for each subtype of lung cancer.It is obvious that in all these three studies, the limited number of lung cancer specimens involved in analysis did not allow for reliable and robust stratified analysis.Furthermore, since the proportions of SqCC-to-AC in the studied group of lung cancer samples in the above mentioned studies were different, a possibility exists that this may have contributed to differences in outcomes concerning the reference gene ranking.Third, in order to synthesize the cDNA strand upon the isolated mRNA, random hexamer primers were used in the two previous studies, while here we used oligo(dT)primers.This also might be a factor possibly contributing to discrepancies between outcomes.From the considerations above a conclusion can be drawn that prior to any gene expression study a proper validation of the reference gene(s) under given specific experimental conditions and setup as well as for given group of individuals involved in study might be required in order to ensure the reliable qPCR data normalization.It may thus imply that any gene previously recommended as a generally suitable reference gene for gene expression studies on a given tissue type may not be considered as really "generally" suitable for that type of tissue and may thus require validation.
Two approaches have been used in order to analyze the reference gene expression stability in our study: NormFinder and geNorm.Although geNorm was criticized recently due to the fact that it basically ranks candidate reference genes based on the similarity of their expression profiles and thus the analysis of candidate genes involved in the same pathway might be misleading (Vandesompele et al., 2002;Andersen et al., 2004), in our study it revealed exactly the same gene ranking as did the model-based NormFinder.Nevertheless, several issues linked to the results obtained by geNorm still remain: 1.Using only geNorm it would be difficult to say whether HPRT1 and RPLP0 were ranked best just because their expression profiles presented the highest correlation coefficient, or because of their really stable expression.NormFinder was shown to be insensitive to correlations in expression profiles, thus it seems reasonable to assume the really stable expression level of these genes.Possibly, the correlation in gene expression plays only a minor role in geNorm analysis, considering the low rank of GAP-DH, a gene significantly correlated with both RPLP0 and HPRT1.2. It is difficult to determine which of the two genes ranked best by geNorm is really the best.This might be a problem, if only one reference gene would be used for a given expression study.3. The involvement of two reference genes for a given expression study according to the geNorm-based ranking would inevitably lead to introduction of a bigger systemic error into qPCR data analysis than the use of best gene pair determined by NormFinder.Being placed on the top two places in reference gene ranking by both programs, HPRT1 and RPLP0 both show the same orientation of their inter-group variation.Since geNorm does not analyze this kind of variation, these genes would be recommended as the best two genes and the qPCR data normalization would not eliminate the inter-group variation in expression levels but introduce a systematic error into the analysis instead (Andersen et al., 2004).The use of HPRT1 and ESD for data normalization as suggested by NormFinder would result in minimization of the inter-group variation of the reference gene expression levels and thus introduce into analysis the least possible systemic error due to opposite orientation of these genes' inter-group variations.4. The assessment of the least number of reference genes for qPCR data normalization is unambiguous since the average pairwise variation V n /V n+1 between two sequential normalization factors (NF n /NF n+1 ) did not fall beneath the 0.15 cut-off value.The authors claim that in a situation when a number of stably expressed genes is subjected to analysis by geNorm, this variation may not reach such low values (dis-cussion with Jo Vandesompele at the geNorm discussion group http://tech.groups.yahoo.com/group/genorm).The six candidate reference genes involved in our study were previously found to present relatively stable expression levels (Liu et al., 2005;Saviozzi et al., 2006), what was also confirmed in our study by the median fold change values not exceeding the value of 3.5 for any gene.This might thus present a rational explanation for the average pairwise variation V n /V n+1 not reaching the cut-off value and the uncertainty concerning the number of reference genes required for qPCR data normalization.In such situations, Vandesompele and coworkers suggest the use of three best ranked genes (Vandesompele et al., 2002;Andersen et al., 2004;Dheda et al., 2004;Radonic et al., 2004).
Concluding, based on our results we recommend the use of HPRT1 as a single reference gene for qPCR data normalization in gene expression studies on healthy and tumor NSCLC tissue.For more accurate normalization, the use of HPRT1 with ESD is recommended.Nevertheless, we emphasize, that the suitability of selected reference genes for a given study must be validated prior to each study.The use of GPADH might not be appropriate for gene expression studies on NSCLC tissue.

Figure 3 .
Figure 3.The geNorm analysis of the candidate reference genes expression stability.A. The stepwise exclusion of least stable genes by calculating the average M value.The average value of expression stability M was determined using geNorm software for each gene and the ranking of candidate reference genes was performed by an iterative algorithm in which the least stable gene with the highest M value was excluded from the next calculation round.Genes on the x-axis are ordered from the left to the right according to their expression stability from the least stable to the most stable one.Average M values: GAPDH: 1.034; POLR2A: 0.959; YAP1: 0.904; ESD: 0.726; RPLP0 & HPRT1: 0.488.B. Determination of the optimal number of reference genes required for qPCR data normalization.The program first calculates the normalization factor NF n from at least two best performing reference genes (n = 2).The variable Vn/n+1 defines the average pairwise variation between two sequential normalization factors (i.e.V3/4 shows the variation of NF 3 in relation to NF 4 ).Values above each bar represent the exact Vn/n+1.

Figure 4 .
Figure 4.The results of the bootstrap analysis.The empirically obtained distributions of ranks shown in percentage of the candidate reference genes as revealed by bootstrapping of the ranking method utilizing the technique of resampling with replacement (size 500).Stable genes show low ranks.Data are given as percentage of bootstrap samples in which a given rank for a given gene was recorded.

Table 1 . Panel of six candidate reference genes analyzed in the study.
a commercially available reference genes; b putative reference genes

Table 2 . Sequences and related information of primers used.
All primer sequences are given in 5'→3' direction; b PCR efficiencies were determined by PCR DataMiner and confirmed by a dilutionseries experiment (see Material and Methods). a

Table 3 . NormFinder analysis results of the candidate reference genes expression stability.
High expression stability is indicated by a low stability value representing the measure of combined inter-and intra-group variation of an individual gene.The best combination of two genes is chosen in order to achieve the minimal systemic error to be introduced into quantitative PCR data analysis during the normalization process.H, healthy non-malignant tissue; T, tumor NSCLC tissue. a

Table 4 . Correlation analysis of the expression levels of six candidate reference genes in healthy normal and tu- mor lung tissue.
Data are presented as correlation coefficients r revealed by the Spearman's rank correlation analysis of the expression levels between each two candidate reference genes.Normal healthy and tumor lung tissue samples were from 21 NSCLC patients (age 64.7 ± 8.0 yrs); *P < 0.01; & P < 0.05.