Introduction
4-Nitroquinoline-1-oxide (4-NQO) treatment induces oral cancer in rats presenting an invasive malignancy in an appropriate anatomic site, which demonstrates considerable histologic and molecular similarity to oral cancers commonly seen in humans (1,2). Therefore, it is a suitable animal model for in vivo evaluation of the efficacy of chemopreventive and therapeutic agents and screening for biomarkers for both tumor progression and treatment. Different strains of rats have been used in this model (3–5). Studies performed in our laboratory (6) demonstrate that invasive oral cancers induced by the administration of 4-NQO develop in 4–6 months after the first exposure to this carcinogenic chemical, and demonstrate highly reproducible incidence and latency patterns. To understand the molecular basis of oral carcinogenesis and identify gene biomarkers and potential chemopreventive targets for future studies, we recently characterized the molecular alteration of the model by performing gene expression microarray using OSCC samples generated from this model. However, in the process of confirming the microarray results by RT-PCR, we found that it would be necessary to identify a reliable reference gene(s) (RG) for normalization.
Quantitative RT-PCR is a commonly used technique for gene expression analyses. However, RT-PCR is greatly affected by multiple factors including RNA integrity, purity, concentration, the presence of inhibitors for RT or PCR in the samples, primers and enzyme efficiencies, DNA contamination, pipetting errors, as well as the choice of RGs for normalization (7). The effects of most of these factors can be minimized by careful operation or using improved approaches. For example, RNA quality and quantity can be assessed using an Agilent bioanalyzer; DNA contamination can be eliminated by treating RNA samples with DNase I. However, RG(s) should be selected with caution. An ideal RG should have a stable basal expression in different tissues, genders, developmental stages, and experimental conditions and should have similar expression levels to the target genes of interest (8). To date, there is no such gene whose expression fulfills these criteria (9). Housekeeping genes have been widely used as RGs, however, the expression levels of housekeeping genes can be affected by both experimental conditions and tissue structure/cellular compositions (10). Thus, the identification of reliable RGs is crucial and should always precede gene expression analyses (8), since gene expression analyses rely on proper normalization to avoid false positive/negative results, which lead to data misinterpretations and even wrong conclusions. Some effort has recently been diverted to the identification of RGs in different tissue types and experimental conditions and different strategies and statistical approaches have been developed (8,11–14). Gene expression comparison between tumors and normal tissues is frequently made for molecular characterization of specific tumors and identification of tumor-specific biomarkers. In the present study, microarray and qRT-PCR approaches were used to screen for potential RGs between tumors and adjacent phenotypically normal tissues. Statistical approaches such as NormFinder (11) and GeNorm (12) that were designed to identify relatively more stable RGs were used to evaluate the stability of potential RG candidates selected based on microarray and RT-PCR data. Hsp90ab1 was successfully established and validated as the best RG in 4-NQO-induced rat OSCC compared to adjacent normal tissues, paving a way for further molecular characterization of this carcinogenesis model.
Materials and methods
Induction of OSCC by 4-NQO in F344 rats and tissue RNA sample preparation
Induction of OSCC by 4-NQO in F344 rats was previously described in detail (15). Briefly, 6–7 week-old male F344 rats received 4-NQO (Sigma-Aldrich, St. Louis, MO, USA) at a concentration of 20 ppm in their drinking water for 10 weeks; after the 10 weeks, the rats received drinking water without added 4-NQO. The study was terminated at 26 weeks after the first day of carcinogen exposure. At necropsy, the tongue from each animal was carefully excised and all gross oral lesions were charted. The tongue was then bisected longitudinally; half of each tongue was fixed in 10% neutral buffered formalin and processed for histopathologic evaluation. Tumor and adjacent normal tissue were carefully separated from the remaining half of each tongue and were snap-frozen in liquid nitrogen and stored at −80°C for use in molecular studies. Histologically confirmed OSCC and adjacent phenotypically normal tissue were selected for molecular analyses. Total RNA was isolated from paired sets of malignant and normal tissues using the RNeasy Mini kit (Qiagen, Germantown, MD, USA) with on-column DNase I digestion according to the manufacturer's instructions.
Microarray analysis
Total RNA isolated from 11 pairs of neoplastic and adjacent normal tongue tissues were individually subjected to microarray analysis using Agilent Rat GE 4x44K v3 arrays (Agilent Technologies, Santa Clara, CA, USA). After RNA isolation, the quality and quantity of total RNA were determined using an Agilent bioanalyzer. First and second strand cDNAs were prepared from the total RNA samples; cRNA target was prepared from the DNA template, verified using the bioanalyzer, fragmented to uniform size, and then hybridized to the microarrays. Slides were washed and scanned using an Agilent G2565 Microarray Scanner. Data were analyzed using Agilent Feature Extraction and GeneSpring GX v7.3.1 software packages. Microarray data for the 11 sample pairs have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GEO), accession GSE51125.
Quantitative RT-PCR
RNA integrity was examined by 1% agarose gel electrophoresis and RNA purity was assessed by spectrophotometer using the A260/A280 ratio. RNA samples with A260/A280<1.8 were eliminated for RT-PCR analysis. RNA was quantitated using Quant-iT RiboGreen RNA assay kit (Life Technologies, Grand Island, NY, USA). RT-PCR analysis was performed as described previously (16) with some modifications. Briefly, 200 ng total RNA was used for RT reaction (in 20 µl reaction volume) for each sample. RT products (cDNA) were diluted by 4-fold with DNase/RNase free water. Real-time PCR was performed with 2 µl diluted RT products in a CFX96 Real-Time PCR detection system using iQ SYBR Green PCR Supermix (both from Bio-Rad, Hercules, CA, USA). Gene-specific primers were designed using Primer 3 software (http://frodo.wi.mit.edu/cgi-bin/primer3/primer3_www.cgi). Primers were designed to span one intron if possible. Primer sequences are provided in Table I. Cq values for each sample were determined by the CFX manager 3.0 software. Since PCR efficiencies were high (>90%) and comparable based on our tests with multiple primer pairs using serial dilutions of pooled cDNA samples, we used the theoretical value of 100% for gene expression calculation (17). Relative gene expression was calculated using the formula 2−(ΔCq), where ΔCq is Cq(sample) − Cq(min). One sample with the highest Cq value [Cq(min)] served as a common control for relative gene expression calculations (8). To ensure the specificity of PCR for each primer pair, melting curve analysis was performed after PCR reaction.
Data analysis
Microarray data sorting and statistical analysis were performed with Microsoft Excel. NormFinder (11) and GeNorm (12) were used to analyze the stability of selected RG candidates in normal and OSCC samples based on our microarray and RT-PCR datasets. NormFinder is a model-based variance estimation approach. It not only takes into consideration the overall intergroup variation (i.e., tumor compared to normal), but also the intragroup variation. NormFinder can analyze expression data obtained through any quantitative method, e.g., real-time RT-PCR and microarray based expression analysis. It ranks the set of candidate normalization genes according to their expression stability in a given sample set and given experimental design.
GeNorm is an algorithm that selects an optimal pair of RGs out of a larger set of candidate genes. It calculates and compares the M-value of all candidate genes, eliminates the gene with the highest M-value, and repeats the process until there are only two genes left. An M-value describes the variation of a gene compared to all other candidate genes. The last pair of candidates remaining is recommended as the optimum pair of reference genes. It is assumed that the candidate genes are not co-regulated (12).
RefFinder (18) was also used to rank RG candidates based on RT-PCR dataset. RefFinder is a web-based interface developed for evaluating and screening RGs from extensive experimental datasets. It integrates the currently available four major computational programs including GeNorm, Normfinder, BestKeeper, and the ΔCt method to compare and rank the tested RG genes. Based on the rankings from each program, it assigns an appropriate weight to an individual gene and calculates the geometric mean of their weights for the overall final ranking.
Results
Induction of OSCC in F344 rats by 4-NQO
Two independent oral cancer induction experiments (28–30 rats per experiment) were performed in F344 rats to generate tissue samples for molecular studies. In these studies, drinking water administration of 4-NQO to F344 rats induced a range of premalignant and invasive malignant lesions in the tongue. The induction of invasive OSCC by 4-NQO was highly reproducible: 83% (25 out of 30) and 75% (21 out of 28) of rats in the two experiments demonstrated invasive oral cancers at six months after the start of carcinogen administration. After histopathologic evaluation, 11 pairs of tissue samples from the first experiment and 16 pairs of tissue samples from the second experiment were selected for microarray and RT-PCR analyses, respectively.
Microarray analyses to screen for potential RG candidates in OSCC compared to normal tissue
Microarray analyses were performed on 11 matched tissue pairs from experiment 1 to compare patterns of gene expression at the mRNA level in OSCC and adjacent phenotypically normal oral tissues. Several approaches were used to select RG candidates from microarray dataset for RT-PCR validation.
Selection of RG candidates based on descriptive statistics
In the normalized microarray data, gene expression levels were expressed in intensity with a range of 0.01–500. Considering that RGs should be easily detectable by RT-PCR, genes expressing at low levels with intensity <5 (intensity cut-off value) were not preferred and therefore removed, which eliminated 90% of the genes in the array list for RG candidates. In addition, genes that were not identified (no name provided or with names such as LOCxxxx or RGDxxxx) were also removed from the array list for RG selection. Then, p-values and fold change between OSCC and normal tissue and CV for the fold change were calculated based on pairwise comparison. Genes with p<0.5 and CV (fold change) >15% were further removed. After these steps, 10 RG candidates were selected and listed in Table II in the order of p-value from largest to smallest. The fold change of these 10 RG candidates between OSCC and normal tissues is close to 1, suggesting that these genes were stable in both OSCC and normal tissues (p>0.5, n=11).
Ranking the selected RGs by NormFinder and GeNorm
Although the best 10 RG candidates were selected from our microarray dataset based on gene expression levels and statistical analysis, it was necessary to rank them based on gene stability to determine the best RG(s). To do so, these RG candidates were further analyzed using NormFinder and GeNorm (Fig. 1) based on microarray data. NormFinder identified Nono as the best RG with a stability value of 0.018 (Fig. 1A), and Nono and Sumo2 as the best combination of two genes with a stability value of 0.013. Dazap2 and Aamp were ranked as the third and fourth most stable RG candidates by NormFinder. GeNorm ranked Dazap2 and Rbm39 as the best pair of RGs with a stability value of 0.1 (Fig. 1B). Nono and Sumo2 were ranked as the third and fourth most stable RG candidates by GeNorm. Apparently, the top RG(s) selected by these programs have not been reported previously.
Expression of commonly used RGs in oral cancer compared to normal tissues
Commonly used RGs such as GAPDH, ACTB are not in the list of the 10 best RG candidates (Table II). Since microarray-based RG selection could be biased, in order to expand the pool of top RG candidates, literature was searched for commonly used RGs in rat and mouse tissue samples and 27 RG candidates were identified (8,19–21) (Table III). Identification of suitable RG(s) for a specific project from commonly used RGs is also a frequently used strategy for RG selection. The expression status of these 27 RG candidates in OSCC compared to normal tissues was analyzed based on microarray dataset. Fifteen out of 27 RG candidates were significantly altered in OSCC compared to normal tissues (p<0.05, n=11) (Table III), suggesting that these 15 genes cannot serve as RGs for OSCC compared to normal tissues, although some of them may serve as RGs for OSCC or normal tongue tissues respectively. Notably, GAPDH was significantly downregulated in OSCC compared to normal tissues; ACTB was significantly upregulated in OSCC. The other 12 RGs were potential RG candidates for OSCC compared to normal tissues. However, Tbp and Tnks were expressed at very low levels and therefore were not suitable RGs for OSCC. Thus, 10 RG candidates were selected for further analysis with NormFinder and GeNorm.
Ranking the commonly used RGs by NormFinder and GeNorm
The selected commonly used 10 RG candidates were ranked by NormFinder and GeNorm (Fig. 2). NormFinder identified Hsp90ab1 as the best RG with a stability value of 0.044 (Fig. 2A), which is higher than the stability value of Nono. However, NormFinder identified Hsp90ab1 and HPRT1 as the best combination of two genes with a stability value of 0.024, which was much lower than the stability value of Hsp90ab1. Consistently, GeNorm also identified Hsp90ab1 and HPRT1 as the best pair of RGs with a stability value (M-value) of 0.18 (Fig. 2B), which is higher than that of the best pair of RGs (Dazap2 and Rbm39) directly selected by GeNorm from the microarray dataset (Fig. 1B).
Determination of top RG candidates and expression of oral cancer biomarker NOS2 and PTGS2 in OSCC compared to normal tissues using top RG candidates for normalization based on microarray data
We selected the top four RG candidates ranked by NormFinder (Nono, Sumo2, Dazap2 and Aamp) and GeNorm (Dazap2, Rbm39, Nono and Sumo2), respectively based on the above microarray dataset analysis (Fig. 1) and top two genes (Hsp90ab1 and HPRT1) selected from commonly used RGs (Fig. 2) for RT-PCR validation. In addition, Glod4 was also selected for RT-PCR validation, since Glod4 was ranked as the least stable gene among the 10 RG candidates by both NormFinder and GeNorm (Fig. 1) and therefore could serve as a reference for comparison. However, while designing primers, we found it difficult to generate very specific primers for Sumo2 and Rbm39 due to the presence of multiple transcript variants/isoforms and these two genes were therefore eliminated for further validation.
We first selected two OSCC biomarkers NOS2 and PGTS2 that were reported to be upregulated in OSCC (15,21–23) to validate the top RG candidates ranked by NormFinder and GeNorm based on microarray dataset. As shown in Table IV, after normalization to these selected RG candidates or combination of two RG candidates, both NOS2 and PTGS2 demonstrated consistent and statistically significant upregulation in OSCC compared to normal tissues, which was also comparable to the microarray results. Consistent with previous analyses (Fig. 1), although Glod4 appeared to be an acceptable RG candidate, it was less stable considering the greater CV and p-values in comparison to other RG candidates (Table IV). These results suggest that these top RG candidates or combination of two RG candidates may be suitable RGs for gene normalization in 4-NQO-induced OSCC in F344 rats. It is still difficult to tell which RG candidate or combination of RG candidates is the best among the selected top RG candidates at this time-point.
Validation of RGs by RT-PCR and determination of the best RG(s)
Ideal RGs should be stably expressed in different sets of samples and validated with a different method. Since these RG candidates were selected based on microarray data and would be used for RT-PCR data normalization, these RG candidates need to be validated by RT-PCR. We therefore used a second discrete set of samples (16 pairs) to validate these top RG candidates. Every effort was made to assure the quality of the RNA samples and RT-PCR. Since 18S rRNA is also traditionally used as an RG in tissue samples, it was included for comparison. Table V shows the analysis of the selected seven RG candidates at the Cq level. Notably, the expression of Dazap2, Aamp and 18S was significantly altered in OSCC compared to normal tissues (p<0.05, n=16) by qRT-PCR and these three genes were excluded for further analysis. Thus, the four remaining RG candidates (Nono, Hsp90ab1, Glod4 and HPRT1) were ranked by NormFinder and GeNorm (Table VI). Hsp90ab1 was the best RG and the best combination of two genes for normalization was Hsp90ab1 and HPRT1 based on NormFinder analysis of the qRT-PCR dataset. GeNorm also consistently identified Hsp90ab1 and HPRT1 as the best pair of RGs. To help further determine the best RG among the four candidates, we used a third program RefFinder to rank these RGs and Hsp90ab1 was still ranked as the best RG (Table VI). Overall, Hsp90ab1 can be considered to be the best RG stably expressed in OSCC and normal tongue tissues in F344 rats treated by 4-NQO. The combination of two genes Hsp90ab1 and HPRT1 can be a better option for normalization with higher stability than that of the single RG-Hsp90ab1 based on NormFinder analysis.
Discussion
With the extensive use of qRT-PCR technique in gene expression analysis, RG identification and selection have become more and more important, since proper normalization of accurate gene expression is crucial for data interpretation and for conclusion to be drawn. In the field of cancer, screening for potential biomarkers by qRT-PCR, is important in determining the RG(s) in normal tissues and cancer at a specific site. In the present study, we identified Hsp90ab1 as a stable RG for normalization of RT-PCR data in 4-NQO-induced OSCC compared to adjacent normal tissues in F344 rats. Similar studies in other cancer types have been previously described (11,17,24,25), however this is the first report of the identification of suitable RG(s) in the 4-NQO-induced oral carcinogenesis model in F344 rats.
Our strategy for the identification of RG(s) in paired OSCC and normal tissue samples is summarized in Fig. 3. This strategy both assures the selection of reliable RG(s) using a combination of two techniques, two sets of RG candidates, two sets of samples and two analytical approaches to avoid some systematic bias, and also minimizes the cost of the study, which is important for research projects, since identification of RG(s) is generally not the major aim in almost any research/clinical projects.
We initially expected that the best RG would be from the RG candidates directly selected from the microarray dataset. However, the best RG after RT-PCR validation turned out to be Hsp90ab1, an RG candidate identified from the literature. This fact suggests that any single approach has limitations, and these limitations should be overcome by using alternate strategies. In addition, although microarray is important for RG screening and selection, RT-PCR validation is indispensable. In fact, some researchers only use the qRT-PCR technique for screening and identification of RGs (12,21); however, if micro-array data are available, this will help to significantly decrease the amount of RT-PCR performed.
It should be pointed out that the criteria established for initial screening is critical for making the RG candidate pool, and further affect the next steps in RG selection. The criteria rely on experimental conditions. In the present study, our aim was to identify RG(s) for OSCC compared to adjacent phenotypically normal tissues; therefore allowing for no statistical difference in RG expression levels between OSCC and the adjacent normal tissues. In this case, p-value became an important parameter for RG selection (p>0.5 was used), which is also used for RG selection in the literature (26). It is possible that some good RG candidates were not identified using the initial screening procedure based on these criteria; however, the RGs selected based on these criteria are acceptable for gene normalization.
No one statistical approach can cover all variables associated with gene expression studies. Thus using more than one statistical approach for RG identification is also an important strategy. The advantage and limitations of different statistical approaches were well discussed in the literature (8). In our studies, NormFinder was preferentially used due to its advantage in considering the variance between subgroups (OSCC compared to normal), GeNorm and RefFinder were also used for comprehensive ranking of the RGs. Hsp90ab1 was selected as the top RG by all these programs in the RT-PCR dataset. Hsp90ab1 has also been previously demonstrated to be stably expressed in human normal and malignant ovarian tissues (27). However, it was regulated by estrogen treatment in mouse uterus (28). Apparently, this RG can be tissue and experiment condition-specific.
Combination of two RGs was more stable than single RGs as identified by NormFinder. We also observed improvement of CV and p-values when the expression of oral cancer biomarker NOS2 was evaluated in OSCC compared to adjacent normal tissues after normalization to the combination of Hsp90ab1 and HPTR1 based on microarray data, suggesting that the combination of two RGs is a choice for gene normalization in tissue samples. However, this improvement was not observed in PTGS2 expression and we believe this was largely due to the intrinsic variation of PTGS2 expression in tissue samples instead of normalization.
More than 50% of the commonly used or traditional RGs were significantly altered in OSCC compared to adjacent normal tissues, demonstrating that RGs must be established for each model and experimental conditions. The biological significance of these alterations is not clear and the structural difference of cell composition in OSCC compared to adjacent normal tissues may account for these alterations. Whether a panel of these housekeeping genes can serve as biomarkers for diagnosis or treatment is of great interest and merits further investigation. In addition, 18S rRNA has been widely used as a stable RG for normalization; however, the RT-PCR data indicate that 18S rRNA expression level was significantly higher in OSCC compared to the adjacent normal tissues; if 18S rRNA is used as an RG for gene expression studies, bias or false results may be generated in this model.
In summary, a simple and relatively economic strategy has been developed for the identification of reliable RGs in cancer compared to normal tissues. The present study identified Hsp90ab1 as the most stable single RG, and Hsp90ab1 plus HPRT1 as the best combination of two genes for normalization in gene expression studies in 4-NQO-induced rat oral carcinogenesis model.