Evaluation and validation of candidate endogenous control genes for real-time quantitative PCR studies of breast cancer

Breast cancer is the most common form of malignancy among women in almost all of Europe and in North America. Each year over one million women worldwide are diagnosed with the disease and it causes over 400,000 deaths annually. Much of the current translational research in this area is based on "transcriptomics", the elucidation of the transcriptional programs underlying disease initiation, promotion and progression, through tumour gene expression profiling. Real-time quantitative PCR (RQ-PCR) [1, 2] is one of the most sensitive and specific quantitation methods for gene expression analysis and is firmly established as a mainstream research tool [3–5]. With the development of high throughput and reliable instrumentation, improved detection chemistries, more efficient protocols and appropriate analysis software, RQ-PCR has become the basis of many breast cancer biomarker studies as well as several novel diagnostic and prognostic assays [6–12]. In addition, RQ-PCR is used to validate microarray expression profiles and quantify genes of interest identified from those analyses.

The most commonly used method to quantify gene expression involves the analysis of target gene expression relative to a control gene. As in other relative gene expression analysis techniques such as Northern blotting and ribonuclease protection assays (RPAs), normalisation of RQ-PCR data is required to control for variation introduced during the steps from RNA extraction to quantitation, especially to control for differences in the quantity and quality of RNA used in reactions [13]. The use of endogenous control (EC) genes, known variously as housekeepers, reference or simply control genes is based on the principle that these genes are expressed constitutively and uniformly in all test samples, so that expression of the target gene can be normalised against them to control for systematic variation in sample handling. Results are then expressed as the ratio of target gene expression relative to the EC gene. In many cases, control genes are inherited from earlier studies using less sensitive forerunner techniques such as Northern blotting and little if any consideration has been paid to validating these genes as controls for specific experiments.

Two of the most commonly used endogenous control genes for breast cancer gene expression studies are glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and β-actin (ACTB) but their reliability in this context has not been demonstrated. In some studies the use of either gene may be inappropriate, as these genes have been implicated in aspects of disease aetiology [14–20]. For example, bisphosphonates; used to inhibit bone resorption in diseases including osteoporosis, Paget's disease and metastatic breast cancer, target GAPDH, decreasing its expression in both breast and prostate cell lines [21], while the actin filament protein family, of which ACTB is a member, may be modulated in malignancy [22] particularly during processes involving reorganisation of the cytoskeleton such as invasion and migration.

The precision of the estimate of change in target gene expression is dependent on the stability of the endogenous control, the variability associated with the target gene and any covariance between the two. Thus the use of non-validated endogenous control genes results in, at best, unreliable data. It is now recognised that a universal, invariably expressed gene is unlikely to exist [23] and may not exist even within individual tissue or cell types. The goal therefore is to identify the most reliable gene or set of genes as endogenous controls for a particular experiment. As a result of this, several groups have developed statistical models and software programs for the analysis of candidate gene stability. The aim of this study was to identify the most stable endogenous control genes from a panel of eleven candidates commonly used as endogenous controls in the context of, but not limited to, breast cancer: GAPDH, TFRC, ABL, PPIA, HPRT1, RPLP0, B2M, GUSB, MRPL19, PUM1, PSMC4, for the quantification of gene expression by relative comparative RQ-PCR in primary breast cancer tissues. The oestrogen receptor alpha (ESR1) transcript, a gene of special significance in breast cancer, was used as a target gene to compare the effect of choice of EC on the estimate of gene quantity.

To identify suitable EC genes for breast cancer gene expression studies in fresh-frozen primary tissue, a panel of 11 genes commonly used as ECs was selected from the literature for analysis of stability: GAPDH, TFRC, ABL, PPIA, HPRT1, RPLP0, B2M, GUSB, MRPL19, PUM1 and PSMC4. Genes were analysed in tumours recovered from patients with benign or malignant breast disease using RQ-PCR. Stability of candidate EC genes was analysed using two statistical analysis tools, geNorm and NormFinder, which employ different statistical models to define the most reliable EC genes for normalisation. The effect of choice of EC gene on target gene expression was analysed using ESR1 as target.

Range of expression of candidate EC genes and ESR1

The candidate ECs displayed a range of C_t values. Mean C_t values per gene and the range of C_t values for each gene are shown in Table 1. Mean C_t values ranged from 19.13 (± 0.21 s.e.m.) for B2M to 26.48 (± 0.15 s.e.m.) for MRPL19.MRPL19 showed the narrowest range followed by PPIA. The genes broadly fell into two categories, those highly expressed with mean C_t values of 19–20 (B2M, RPLP0, GAPDH, PPIA) and moderate abundance genes with mean C_t values of 23–26 (PSMC4, ABL, GUSB, TFRC, PUM1, HPRT1 and MRPL19). The target gene ESR1 showed the broadest range of C_t values (10.35) from 17.10–27.45.

Gene symbol	Ct Range	CtMin.	CtMax.	Mean Ct ± s.e.m.
MRPL19	4.06	24.79	28.84	26.48 ± 0.15
PPIA	4.4	18.01	22.41	19.74 ± 0.18
B2M	4.67	17.22	21.89	19.13 ± 0.21
PUM1	5.06	22.38	27.44	24.78 ± 0.21
GUS	5.19	21.37	26.55	23.94 ± 0.20
PSMC4	5.69	20.84	26.52	23.32 ± 0.22
ABL	5.89	21.27	27.15	23.75 ± 0.30
TFRC	6.03	21.4	27.42	24.40 ± 0.29
GAPDH	6.88	16.79	23.67	19.71 ± 0.25
HPRT1	7.05	21.57	28.61	24.84 ± 0.27
RPLPO	8.8	16.1	24.9	19.46 ± 0.33
ESR1	10.35	17.10	27.45	21.96 ± 0.34

With conversion of C_t values to relative quantity values (Q_Rel.), there was no difference in candidate EC gene quantities between benign and malignant tissues (P > 0.05; Fig. 1a). There was however, a significant difference in variance between genes (P = 0.001; Fig. 1b) indicative of differing stabilities of the candidates.

Analysis of EC gene stability

The significant difference in EC variability demonstrated the necessity to validate their use in this context. Expression stability was analysed using the two softwares geNorm [24, 25] and NormFinder [26].

GeNorm uses a pair-wise comparison-based model to select from a panel of candidate EC genes, the gene-pair showing least variation in expression ratio across samples. It calculates a measure of gene stability (M) of each gene based on the average pairwise variation between all tested genes. Genes with the lowest M values are those demonstrating most stable expression. Table 2 shows the M values for all tested genes. Eight of the eleven genes analysed; PPIA, MRPL19, GAPDH, PUM1, B2M, HPRT1, PSMC4 and TFRC, showed M values less than the geNorm default threshold of 1.5, while the three remaining genes; RPLP0, GUSB and ABL, showed M values greater than that threshold. In a stepwise progression, geNorm excludes the least stable gene, recalculating M for the remaining genes, resulting in the characterisation of the stability of each gene on a ranked scale and ultimately the identification of the two most stably expressed genes (Fig. 2a). As shown,ABL and GUSB were the first and second genes respectively excluded from the analysis on the basis of instability and MRPL19 and PPIA were identified as most stable gene-pair.

	NormFinder		geNorm
Rank	Gene	Stabilityb	Gene	Stability
1	MRPL19	0.105	PPIA	1.116
2	PPIA	0.119	MRPL19	1.134
3	B2M	0.169	GAPDH	1.240
4	GAPDH	0.176	PUM1	1.275
5	PUM1	0.222	B2M	1.305
6	HPRT1	0.233	HPRT1	1.357
7	TFRC	0.236	PSMC4	1.419
8	PSMC4	0.247	TFRC	1.437
9	GUSB	0.268	RPLPO	1.549
10	RPLPO	0.295	GUSB	1.568
11	ABL	0.413	ABL	2.070

^bHigh expression stability is indicated by a low stability value. For geNorm, stability is based on an estimate of the pairwise variation (M) relative expression while for NormFinder, stability is based on combined inter- and intra-group variation.

GeNorm also calculates a normalisation factor (NF) required to determine the optimal number of EC genes required for accurate normalisation. This factor is calculated using the variable V as the pairwise variation (V_n/V_{n + 1}) between two sequential NFs (NF_n and NF_{n + 1}). To meet the recommended cut off V-value of 0.15, the point at which it is unnecessary to include additional genes in a normalisation strategy [24], the programme indicated the use of 5 of the six most stable genes i.e., MRPL19, PPIA, GAPDH, PUM1 and B2M (Fig. 2b). However, there was no significant effect on relative quantity of ESR1 expression using the 5 gene panel of MRPL19, PPIA, B2M GAPDH and PUM1, compared to the two-gene combination of MRPL19 and PPIA (P > 0.05).

Stability of gene expression was also analysed using NormFinder [26]. This programme uses a combined estimate of the intra- and inter-group variation to determine the most stably expressed candidate EC gene and gene-pair. Table 2 shows the ranking of the candidates. As for geNorm, NormFinder identified MRPL19 and PPIA as the most stable pair of genes and MRLP19 as the single most stable gene.

Effect of EC on ESR1 relative gene expression

There was no effect of choice of EC on the relative quantity of ESR1 (P > 0.05) probably due to the large variation associated with the expression of ESR1 compared to any of the EC genes as shown in Fig. 3. However, there was a significant effect of EC gene on the estimate of the error associated with relative gene expression (P < 0.05). The magnitude of the error was significantly reduced using MRL19 and PPIA as a combined endogenous control compared to the use of ABL, B2M, GAPDH, GUSB, HPRT1, PSMC4, PPIA, PUM1, RPLP0 and TFRC. While there was no difference in error using MRPL19 and PPIA in combination compared to MRPL19 on its own (P > 0.05), there was a trend towards a reduction in the estimate of the error using the combined normalising factor. There was a significant increase in relative quantity of ESR1 in malignant versus benign breast tissues compared to normal breast tissue (P < 0.01; Fig. 4).

To our knowledge this is the first systematic evaluation of the reliability of a large number of genes used as endogenous controls for RQ-PCR analysis in breast cancer studies. The literature cites just two articles in relation to the evaluation of EC genes for breast cancer [27, 28]. However, one study characterised ECs not in primary breast tissues but in cell lines [27] and the other, while employing primary breast tissues, compared just two commonly used EC genes with other genes selected from a microarray dataset of breast cancer tissues and cell lines, as well as cell lines of different origins [28].

There was no effect of tissue type on EC expression, validating comparison of their stability. This is an essential but often overlooked precursor analysis when using geNorm and NormFinder [42] since these methodologies assume the candidates are not differentially expressed between experimental groups. There was however a significant difference in variance between candidates (P = 0.001; Fig. 1), with genes such as RPLP0, TRFC, HPRT1 and GAPDH showing greater variance than others e.g., MRPL19 and PPIA. Since the resolution of RQ-PCR is defined by the variance associated with the EC [13] these results emphasise the necessity to evaluate and validate EC genes.

A single universal EC is unlikely to exist [43] and since the function of most genes is largely unknown it is impossible to predict their expression under different experimental conditions. The use of more than one EC hedges the bet and increases the accuracy of quantitation compared to the use of a single EC [13, 24, 26, 36, 44]. Studies show substantial errors, up to 6.5-fold, in expression quantitation using single as opposed to multiple EC genes [24]. In this study, stability of expression was analysed using two distinct statistical models, a pairwise comparison model, geNorm, and an ANOVA-based model, NormFinder. The geNorm applet selects from a panel of genes, the pair showing least variation in expression ratio across samples and estimates the minimum number of genes required for optimal normalisation. NormFinder estimates stability values for ECs considering combined intra- and inter-group variation and identifies the most stable gene and gene-pair, where the stability of the pair exceeds that of the single gene. Despite their differences both models identified MRPL19 and PPIA as the most reliable ECs while RPLP0, GUSB and ABL were least reliable. This result reflects those of the equality of variance analysis and, broadly, ranking by range of C_t values (Table 2). However, the ranking of genes by C_t range and by the model-based methods differed for some genes e.g., GAPDH, demonstrating the necessity to scale and correct raw C_ts for amplification efficiency before analysis. GeNorm indicated that optimal normalisation could be achieved using the five most stable genes but there was no difference in ESR1 gene expression using this approach as opposed to the two-gene combination. While it is not known whether this would hold for other less variable target genes, cost and sample availability are limiting factors for most studies so the two-gene combination may be more practical for most applications.

The effect of using less stable ECs was assessed using ESR1 as a target gene. Due to the high variability of the expression of this gene (Fig. 4) there was no effect of EC on quantitation. Without further discussion of the relevance of the differential expression of ESR1 in benign and malignant tissues; a gene whose role in breast cancer is widely appreciated, it is likely that had the target gene shown a more discreet change in gene expression, an effect of EC on quantitation would have been apparent. There was however a clear tendency for ESR1 expression to change depending on EC (Fig. 4). In the benign samples ESR1 could be made to appear up- or down-regulated depending on EC, while its expression in the malignant samples could be numerically altered by one order of magnitude – artifactual results due simply to the choice of EC. Clearly this type of error is unacceptable, especially in the analysis of markers for potential clinical application.

Furthermore, there was a significant effect of EC on the magnitude of the error associated with the estimate of ESR1 expression. The use of the MRPL19, PPIA combination minimised the error compared to all other ECs with the exception of MRPL19. Apart from the fact that the use of single EC genes can compromise data as already described, the use of either gene alone is not recommended for two further reasons. Firstly, geNorm identified them based on a pairwise comparison model so their individual use is inconsistent with that analysis method. As shown in Fig. 4, ESR1 expression can be made numerically increase or decrease in the benign groups if these genes are used singly. Secondly, stability, as assessed by NormFinder, increased roughly 30 percentage points, from 0.105 using MRPL19 alone to 0.072 using it in combination with PPIA indicating improved reliability of the two-gene combination [26].

One previous study analysed EC stability in primary breast tissue [28]. That study compared two traditionally used ECs, GAPDH and ACTB and four genes identified in microarray studies [45, 46]: MRPL19, PUM1, PSMC4 and SF3A1. The authors recommended MRPL19 be used with PSMC4 and PUM1. However, in the present study, PPIA, which was not assessed by Szabo and colleagues, showed greater stability than PSMC4 and PUM1. Of the genes analysed in this study, PPIA ranked first and second by Normfinder and geNorm respectively compared to fifth and seventh for PSMC4 and eighth and fourth for PUM1. The reason(s) PPIA did mot emerge as a candidate in Szabo and colleagues' microarray study is unclear. As part of their selection procedure, data was filtered to remove near background signals from low abundance genes yet in this study PPIA showed the third highest mean C_t value. Cohort-specific effects are also unlikely since the tissues are broadly similar in terms of their clinical and histological parameters. Apart from the obvious differences in the quantitative capacities of microarray and RQ-PCR technology, a possible explanation is that the cohort from which these authors selected their candidates was not breast cancer-specific and included tissue from metastatic breast cancer lymph nodes and cell lines from an assortment of origins including dermal, leukemic, umbilical and melanoma samples [46].

This study also analysed associations between genes. An assumption of the geNorm model is that candidate ECs are not co-regulated yet the analysis of such genes by that method would lead to an erroneous choice of best gene-pair. To our knowledge the EC genes evaluated here are functionally independent as shown in Table 4. In addition, covariance of target and EC is clearly unacceptable. However, regression analysis showed significant linear relationships between ESR1 and three of the candidates: TFRC, HPRT1 and ABL. Approximately two thirds of breast tumours are oestrogen-dependent and the number of genes whose expression is known to be, or likely to be mediated through the receptor is sizeable [47]. This makes analysis of associations between oestrogen-responsive target and control genes of particular importance in breast cancer studies. Although there is little evidence that these genes are regulated by oestrogen, the results suggest that aside from their poor stability it would be wise not to use them as ECs in the analysis of oestrogen-responsive breast cancer.

Despite the clear increase in accuracy afforded by the use of more than one validated EC gene, a recent survey of working practices indicated that over half of those polled continue to use one reference gene and that two thirds of these do not validate that gene [48]. Obviously levels of awareness and/or willingness to address this issue must be improved.

The current emphasis on personalised cancer treatment has resulted in the development of prognostic and predictive multi-gene RQ-PCR assays. However, with such developments comes the demand for greater accuracy and resolution of gene expression quantitation. In this study we have validated two genes, MRPL19 and PPIA as EC candidates for RQ-PCR analysis of primary breast tissue using two different statistical models and demonstrate that of the genes in current use in this field, the above gene combination offers increased accuracy and resolution in the relative quantitation of gene expression data. The genes identified should be of use in a broad spectrum of translational research and diagnostic applications in breast cancer.