Identification of Suitable Reference Genes for Real-Time PCR Analysis of Statin-Treated Human Umbilical Vein Endothelial Cells

Proper data normalization in quantitative real-time reverse-transcription polymerase chain reaction (RT-qPCR) is of critical importance for reliable mRNA expression analysis. Due to a diversity in putative reference genes expression stability in different in vitro models, a validation of an internal control gene should be made for each particular tissue or cell type and every specific experimental design. A few approaches have been proposed for reference gene selection, including pair-wise comparison approach and model-based approach. In this article we have assessed the expression stability of eight putative reference genes: ACTB, B2M, GADD45A, GAPDH, HPRT1, PES1, PSMC4, YWHAZ, in human umbilical vein endothelial cells (HUVEC) treated with different statins and with TNF-α. The analysis was performed with three reference gene validation programs: geNorm, NormFinder and BestKeeper. We have shown that hypoxanthine phosphoribosyltransferase 1 gene (HPRT1) and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta polypeptide gene (YWHAZ) are the most stably expressed genes among the analyzed ones. Furthermore, our results show that β-actin gene (ACTB) is downregulated by statins and thus should not be used as a normalizing gene in a discussed experimental setup. A ranking of candidate reference genes stability values is provided and might serve as a valuable guide for future gene expression studies in endothelial cells. This is the first report on reference gene selection for RT-qPCR applications in statin-treated HUVEC model.


Introduction
Quantitative real-time reverse-transcription polymerase chain reaction (RT-qPCR) has become one of the most popular techniques of quantifying mRNA levels. The method is relatively easy and precise, however, when used in an inappropriate way, it can lead to considerable misinterpretation of results [1][2][3][4]. It is therefore of critical importance to perform proper data normalization which enables to control differences between samples that may arise at different stages throughout the procedure. However, although the qPCR has become very popular, data normalization still remains a problem. There are several strategies which can be applied to normalize qPCR results [3], but the most common one is the use of reference genes as an internal standard. Although recent studies clearly show the importance of a proper choice of a reference gene [2,5,6], still many currently published reports present RT-qPCR results that miss information on a reference gene selection. In addition, researchers often routinely use the most classical reference genes, such as genes coding for glyceraldehyde 3-phosphate dehydrogenase (GAPDH) or b-actin (ACTB), convinced that these are the universal reference genes and unaware that they can be highly regulated [7][8][9].
The ideal reference genes should be expressed at the same level in all cells and under all experimental conditions. It has been, however, well documented, that most of them undergo significant regulation and thus cannot be considered as a proper reference [10]. Despite these limitations, the use of reference genes as internal controls remains the most common method used to normalize cellular mRNA content in analyzed samples [11,12]. A fortiori, the use of this method should be preceded with rigorous reference genes validation to avoid an improper gene choice. Therefore, for each particular tissue or cell type and specific experimental designs, a thorough search is needed to ensure that no significant change in a reference gene expression occurs [13]. Unfortunately, despite its importance, this experimental step is often neglected.
The aim of this study was to identify the most stable reference genes for human umbilical vein endothelial cells (HUVEC) treated with different 3-hydroxy-3-methylglutaryl-coenzyme A (HMG-CoA) reductase inhibitors, known as statins. HMG-CoA reductase inhibitors are among the most frequently prescribed drugs for the prevention and treatment of cardiovascular diseases. Besides lowering the plasma cholesterol concentration, they exert pleiotropic effects that are independent of their cholesterol-lowering properties [14,15], which include improved endothelial functions and decreased vascular inflammation [14]. Although HUVEC are often used as an in vitro model to determine mechanisms of statins effects on endothelial cells and expression of various genes is analyzed using RT-qPCR method, very little is known about the stability of potential reference genes used in such studies [16]. We have assessed the expression stability of eight putative reference genes in HUVEC treated with six different statins: lovastatin, atorvastatin, fluvastatin, simvastatin, pravastatin and cerivastatin and additionally stimulated with a proinflammatory cytokine, tumor necrosis factor alpha (TNF-a). The reference genes examined were: b-actin (ACTB), glyceraldehyde 3-phosphate dehydrogenase (GAPDH), b-2-microglobulin (B2M), growth arrest and DNA-damage-inducible protein alpha (GADD45A), hypoxanthine phosphoribosyltransferase 1 (HPRT1), pescadillo homolog 1 containing BRCT domain (zebrafish) (PES1), proteasome (prosome, macropain) 26S subunit, ATPase 4 (PSMC4) and tyrosine 3monooxygenase/tryptophan 5-monooxygenase activation protein, zeta polypeptide (YWHAZ). We believe that the presented data related to a reference gene selection will be useful not only for the RT-qPCR analyses of statin-treated HUVEC, but also for other studies with human primary cells.

Experimental Design
In the present study eight candidate reference genes (ACTB, B2M, GADD45A, GAPDH, HPRT1, PES1, PSMC4, YWHAZ) ( Table 1) were evaluated in HUVEC which underwent statin or combined statin and TNF-a treatments. All the necessary controls were included. All HMG-CoA reductase inhibitors, i.e. lovastatin, atorvastatin, fluvastatin, simvastatin, pravastatin and cerivastatin, were used at concentrations which did not induce any cytotoxicity (data not shown). All the selected candidate reference genes belong to different functional classes to minimize the chance of their coregulation.
Data were analyzed in three reference gene validation programs: geNorm [17], NormFinder [18] and BestKeeper [13], and the results were used to rank the candidate reference genes from the most to least stable. Based on the rankings obtained from each program, which assigned appropriate weights to every individual gene, the geometric mean of their weights was calculated for the overall final rankings. Candidate reference genes were analyzed thrice using samples from three different cell donors. 20 cDNAs containing statin and statin-and-TNF-a samples were obtained from each donor. The obtained data were analyzed for each donor either separately for statin-treated cells and statin-and-TNF-atreated cells or in the pooled analysis of all 20 samples.

Cell Culture and Treatment
HUVEC were purchased from Invitrogen Life Technologies. Three populations from different cell donors were used. Cells were grown in EBM-2 basal medium supplemented with the EGM-2 SingleQuots kit (Lonza, USA). For all experiments HUVEC at passage four were used.
Cells were grown to confluence and then treated with lovastatin at final concentrations 1 mM and 2 mM, atorvastatin 1 mM and 2 mM, fluvastatin 1 mM and 2 mM, simvastatin 1 mM, cerivastatin 0.1 mM or pravastatin 1 mM for 48 hours in full EGM-2 medium. Untreated cells were used as controls. The same experimental setup was repeated and followed with TNF-a treatment (10 ng/ mL, 1 hour).

RNA Isolation and cDNA Synthesis
Total RNA was extracted from HUVEC and purified using the NucleoSpin RNA II Kit (Marcherey-Nagel, Germany) according to the manufacturer's instruction. Cell lysis was performed in a RNases inactivating buffer provided by the manufacturer. Until RNA extraction samples were stored at 270uC. Purified RNA was reverse transcribed immediately after extraction.
cDNA was synthesized using a High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, USA) following manufacturer's instructions. The reaction was set with 6 mg of total RNA in a total volume of 60 mL containing: random primers, 4 mM deoxyribonucleotide triphosphates (dNTP's), 2.5 U/mL Multi-Scribe Reverse Transcriptase and RT buffer. RNase inhibitor (1 U/mL) was used for each reverse transcription PCR reaction. Cycle parameters were set to 10 minutes at 25uC, 120 minutes at 37uC and 5 minutes at 85uC. cDNA was stored at 220uC until further use.

Quantitative Real-time Reverse-transcription PCR
The following eight putative reference genes were selected for analysis: ACTB, B2M, GADD45A, GAPDH, HPRT1, PES1, PSMC4 and YWHAZ. The selected genes belong to different functional classes, which reduces the chance of co-regulation.
All primers and probes were purchased from Applied Biosystems, USA ( Table 2). The real-time qPCR reactions were performed using TaqManH Gene Expression Assays (FAM TM dye-labeled MGB probes) and TaqManH Gene Expression Master  PCR conditions were as follows: 50uC for 2 minutes, 95uC for 10 minutes followed by 40 cycles of 95uC for 15 seconds (denaturation step) and 60uC for 60 seconds (annealing and extension step) during which fluorescence was measured. Data expression levels were recorded as quantification cycles (C q ). Data was acquired using the 7500 Software (Applied Biosystems, USA). The mean C q values of the triplicate reactions were used in further analysis.
Calculations of the mean C q and fold change values were performed by means of DataAssist v3.0 Software (Applied Biosystems, USA).

PCR Efficiency
A 2-fold dilution series was prepared from pooled cDNA samples. PCR reactions were performed as described above in triplicates. The PCR efficiencies (E) and correlation coefficients (R 2 ) for each primer pair were calculated using the formula: E (%) = (10 (-1/slope) -1)N100. Slopes were determined from a standard curve obtained when a logarithm of the initial template concentration was plotted on the x axis and C q on the y axis [11]. The efficiencies for all primer pairs are listed in Table 2.

Data Analysis
GeNorm analysis. GeNorm is a software for Microsoft Excel which provides a measure of gene expression stability [17]. It ranks the genes basing on the internal control gene stability parameter M. M is the mean pair-wise variation between individual gene and the other putative reference genes tested. Stepwise exclusion of the gene with the highest M value and recalculation allows ranking of the tested genes according to their expression stability. Lower M values represent higher expression stabilities. Any gene with M .1.5 is suggested to be considered unreliable as a stable reference gene [17]. GeNorm authors suggest the use of 10 samples and 8 reference genes for validation procedure.
Normfinder analysis. NormFinder is another Excel-based statistical algorithm that computes expression stability values to range candidate reference genes [18]. A high stability value represents a high gene expression variance. In addition, the program allows for comparison of inter-and intra-group variation of gene stability. NormFinder authors suggest using at least 8 samples per group and minimum 3 candidate genes, but recommend 5-10 genes [18].
BestKeeper analysis. BestKeeper analyses C q values to evaluate the expression variability of the reference genes. The key factor in the analysis is the standard deviation (SD) which represents the stability of the gene, and the lower SD value, the better stability. Any studied gene with the SD .1 is suggested to be considered unreliable [13].
Afterwards, the program performs a comparative analysis based on pair-wise correlation coefficient (r) between each gene and the BestKeeper Index (BI), which is the geometric mean of C q values of candidate reference genes.
The genes with SD .1 are eliminated from further analysis and the remaining genes are ranked according to their coefficient of correlation (r).
Final ranking. All analyzed genes were ranked by all three programs and appropriate weights to every gene were assigned. For the overall final rankings the geometric means of the obtained weights were calculated.

Results
For each primer pair PCR efficiencies were calculated from the slope of the standard curve. The obtained efficiencies varied from 90.5% to 101.1% (Table 2).

First Donor Assay
The first assay was performed with a set of eight primer pairs (ACTB, B2M, GADD45A, GAPDH, HPRT1, PES1, PSMC4, YWHAZ) for 20 cDNA samples, i.e. statin and combined statinand-TNF-a-treated cells from the first donor. C q values were used for further analysis. Three reference gene validation programs were used: geNorm, NormFinder and BestKeeper. Mean C q values were input into BestKeeper. For geNorm and NormFinder C q values were transformed into relative quantification data using the equation 2 (-DCq) . DC q is the difference between data point of interest and the highest data point in the data set. Therefore all data is relative to the sample showing the lowest level of gene expression. For pooled analysis of 20 cDNA samples, including statin-and combined statin-and-TNF-a-treated samples, GeNorm ranked analyzed genes basing on their stability value (M). The most stably expressed genes, with the lowest M-value, were HPRT1 and YWHAZ (Table 3, part A). NormFinder also ranges genes depending on a stability value but using different algorithm. The lowest stability value was for B2M and the second lowest value was for GAPDH (Table 3, part A). BestKeeper produced descriptive statistics (data not shown). For all the analyzed genes, except for ACTB, SD values were below 1 ( Table 3, part A), which suggests that ACTB should be excluded from further analysis and the rest of analyzed genes could be potentially used as reference genes. The ranking of putative reference genes based on coefficient of correlation values (r) is shown in Table 3, part A. The best correlated genes were YWHAZ followed by HPRT1. Next, geometric means of the weights from all three rankings for every individual gene were calculated and are presented in Table 3, part A. In the final ranking YWHAZ is positioned first and is followed by HPRT1.
The statin-and combined statin-and-TNF-a-treated samples were also analyzed separately. For statin-treated samples the most stably expressed genes according to geNorm analysis were HPRT1 and B2M. NormFinder ranked HPRT1 at the first and YWHAZ at the second position. According to BestKeeper the best correlated genes were YWHAZ and HPRT1. In the final ranking HPRT1 was positioned first and YWHAZ second. The results are presented in Table 4, part A.
For combined statin-and-TNF-a-treated samples geNorm analysis indicated HPRT1 and B2M as the best reference genes. NormFinder ranked YWHAZ as the best reference gene followed by PSMC4. The best correlated genes according to BestKeeper were HPRT1 followed by GAPDH. Final ranking indicated HPRT1 as the best reference gene and YWHAZ as the second best reference gene (data not shown).

Second Donor Assay
The second donor assay was also performed with 20 cDNA samples and the eight primer sets evaluated in this study. The methods of analysis were the same as previously. For pooled analysis of 20 cDNA samples, including statin-and combined statin-and-TNF-a-treated samples, the most stably expressed genes according to geNorm are HPRT1 and YWHAZ. NormFinder indicated B2M and YWHAZ as the best and the second best reference genes respectively. BestKeeper analysis ranked YWHAZ as the best reference gene followed by PSMC4. In the final ranking, similarly to the first cell donor, YWHAZ and HPRT1 were ranked as the best and the second best reference genes respectively. The results are presented in Table 3, part B. Again, samples obtained from statin-treated cells were analyzed separately from combined statin-and-TNF-a-treated samples. All respective controls were included in the analyses. The most stably expressed genes for statin-treated cells according to geNorm analysis were HPRT1 and YWHAZ. NormFinder analysis suggested YWHAZ as the best and PSMC4 as the second best reference genes. HPRT1 was ranked at the third position. BestKeeper analysis ranked YWHAZ at the first and HPRT1 at the second position. For all analyzed genes SD values were below 1. Final ranking ranked YWHAZ at the first position followed by HPRT1. The results are presented in Table 4, part B.
For samples treated with statins and TNF-a geNorm ranked B2M and YWHAZ as the two best reference genes. NormFinder positioned PSMC4 and YWHAZ at the first and second position respectively. According to BestKeeper analysis the best reference gene was HPRT1 followed by B2M. In the final ranking B2M and PSMC4 were ranked at the first two positions. YWHAZ and HPRT1 were ranked at the third and fourth position respectively (data not shown).

Third Donor Assay
Similarly to the first and second donor, the third donor assay was performed with eight primer sets and 20 cDNA. The methods of analysis were the same as previously. For pooled analysis of 20 cDNA samples, including statin-and combined statin-and-TNF-atreated samples, the best reference genes indicated by geNorm analysis were B2M and YWHAZ. NormFinder ranked PSMC4 at the first and HPRT1 at the second position. According to BestKeeper the best correlated reference genes were PSMC4 and HPRT1. In the final ranking PSMC4 was ranked first and HPRT1 second. YWHAZ was positioned fourth. The results are presented in Table 3, part C. In a manner analogous to the first and the second donor assays, results obtained for statin-treated samples were analyzed separately from combined statin-and-TNF-a-treated samples. All respective controls were included in the analyses. GeNorm analysis indicated YWHAZ and PSMC4 as the two best reference genes for RT-qPCR studies with statin-treated HUVEC. NormFinder and BestKeeper ranked HPRT1 as the best and B2M as the second best reference gene. In final ranking HPRT1 was ranked first and YWHAZ second. The results are presented in Table 4, part C.
For combined statin-and-TNF-a-treated samples geNorm ranked HPRT1 and GAPDH as the two best reference genes. According to NormFinder and BestKeeper GAPDH and PSMC4 were the best and second best reference genes respectively. At the final ranking GAPDH was positioned highest and was followed by HPRT1. YWHAZ was ranked fourth (data not shown).

Discussion
RT-qPCR has become a gold standard for quantifying mRNA. As this method of analyzing gene expression is highly specific and relatively easy, it has reached a great popularity. However, the data normalization still remains an issue. The most common method of normalizing qPCR results is the use of reference genes and the strategy is based on the assumption that they are stably expressed. However, many studies have demonstrated that the stability of each reference gene needs to be verified individually under specific experimental conditions [7][8][9]. It is also recommended that a series of genes are tested for stability and more than one is used for normalization in the final experimental setup. When a gene of interest is not compared to appropriately validated, stably expressed reference genes, misinterpretation of results may occur. Constantly growing evidence indicates that there is no single reference gene that can be used for different experiments, but hopefully with the growing number of experimental data and reports, such as this one, a group of putative reference genes for certain specific experimental setups could be recommended for future studies [5,6,19,20].
In this study we have shown the variability in the expression stability of eight putative reference genes (ACTB, B2M, GADD45A, GAPDH, HPRT1, PES1, PSMC4, YWHAZ) in statin-treated HUVEC when compared in three Excel-based programs: geNorm, NormFinder and BestKeeper. As for the validation of a reference gene only limited number of samples from all to be analyzed are usually used, we decided to verify the impact of the samples' selection on the obtained results. For this purpose we performed three assays with three sample sets collected from different cell donors. We have also verified weather the reference genes selected for statintreated cells may also be used for RT-qPCR analysis of cells additionally stimulated with TNF-a. Therefore, for each cell donor three sets of analyses were performed: 1) for statin-treated cells, 2) for combined statin-and-TNF-a-treated cells, and 3) pooled analysis for all the samples.
Although the results obtained for every donor differ to some extent, certain putative reference genes (i.e. HPRT1, YWHAZ and B2M) are ranked high in most of the analyses, while the other (i.e. ACTB, GADD45A and PES1) are predominantly ranked low.
In all the analyses of statin-treated cells, ACTB, one of the most popular reference genes, has been ranked at the last position (Table 4) suggesting that statins affect its expression. A more detailed evaluation of the obtained results has indicated that the expression of ACTB is down-regulated by all statins (data not shown), what excludes it from a group of potential reference genes for the presented type of experiments.
GAPDH, another reference gene often used to normalize RT-qPCR data without any validation, has been ranked low indicating that it is not suitable for this research model. These results show that the validation of a reference gene for normalizing RT-qPCR data is crucial and using popular reference genes, such as ACTB or GAPDH, without any validation may lead to false results.
Summarizing, HPRT1 and YWHAZ were ranked high in all the analyses which makes these genes the best choice for normalizing gene expression in statin-treated HUVEC. As it is commonly suggested to use more than one reference gene for normalizing data in qPCR studies [3], we recommend these two genes as reference genes in the presented experimental setup.
The differences in the results obtained from the first, second and third donor assays may reflect a normal genetic diversity of human population. HUVEC are primary cells isolated from human umbilical vein and in this study each pool has originated from a different donor. This might be the reason for some differences in the genes regulation in the presented experimental setup leading to the disparity of the rankings. Therefore, our study shows that the validation of reference genes for experiments based on heterogeneous cell cultures requires more samples than recommended by geNorm and NormFinder authors.
The disparities between the programs output in each analysis are a result of different methodologies used in the calculation of gene stabilities (e.g. model-based approach and pair-wise variation evaluation) and overall inherent variability of the genes examined. The pair-wise comparison approach (geNorm, BestKeeper) selects the most suitable reference gene on the basis of the variation of expression ratios between candidate reference genes expression across the sample set. It is based on the assumption that the ratio between two putative reference genes is constant across samples independently to RNA amount analyzed per sample. The variation of this ratio for two candidate reference genes across samples (pair-wise variation) is a measure of gene stability. However, geNorm and BestKeeper algorithm analyses are based on the assumption that none of the genes analyzed in the study is coregulated.
The co-regulation of candidate reference genes does not significantly affect the model-based approach (NormFinder). Nevertheless, this type of analysis can be sensitive to sampling errors and outliners. For that reason the use of more than one type of algorithm for the validation of reference genes is suggested. The comparison of the reference genes rankings obtained from more than one program will give more reliable results.
This report should draw particular attention to a proper experimental planning. In the optimal setup reference gene validation should be carried out for every experiment and every pool of samples, but as the process is time and money consuming this recommendation seems difficult to achieve. Nevertheless, our study clearly shows that the more samples used for the validation of the reference gene the better. The same rule applies to choosing reference genes, however, the possibility that some of the selected genes are co-regulated and thus may falsify the results should be considered.
In conclusion, in this study we have shown that among the analyzed genes, HPRT1 and YWHAZ are the most suitable reference genes for the expression studies in HUVEC treated with statins and additionally stimulated with TNF-a. Moreover, our results clearly show that ACTB should not be used as a normalizing gene in a discussed experimental setup. These data may also be useful when validating reference genes for other studies with HUVEC. Our observations confirm that the proper selection of a reference gene is crucial for reliable data analyzing.