Normalization with genes encoding ribosomal proteins but not GAPDH provides an accurate quantification of gene expressions in neuronal differentiation of PC12 cells

Selection of candidate reference genes from microarray data

Real-time PCR validation of novel candidate reference genes

As one of the most extensively studied models for neuronal differentiation, PC12 cells respond to a broad spectrum of pharmacological agents, which trigger a myriad of intracellular signaling pathways leading to neuronal differentiation. In order to verify the general utility of the 20 selected putative reference genes (Table 1) in a broader range of experimental conditions, we differentiated PC12 cells with other stimuli (Forskolin [33], KCl [34] and ROCK inhibitor Y27632 [35]) in addition to NGF and GDNF. GDNF was applied to PC12 cells stably expressing GDNF Family Receptor alpha 1a (GFRa1a) and co-receptor RET (either RET9 or RET51 isoforms), which are not endogenously expressed at detectable levels in PC12 cells (data not shown). The percentage of PC12 cells differentiated by the five chemical stimuli was quantified (Figure 1A) and the axon-like features of the extended neurite were confirmed by immunocytochemical analysis with anti-Neurofilament-200 antibody (Figure 1B). The extent of neurite outgrowth was highly dependent on the stimuli used. NGF and GDNF stimulation induced longer neurite outgrowths than Forskolin, KCl or Y27632. Total RNA was collected at 0.5 h, 6 h, 24 h and 72 h from control and treated cells for each stimuli, with biological triplicates, that totaled 120 samples.

We analyzed the expression levels of the aforementioned twenty candidate reference genes, two most commonly used HKGs (GAPDH and ACTB), and three genes which are well known to be regulated by NGF using RT-qPCR (Figure 2). The expression levels of the twenty candidate reference genes and the two HKGs span three orders of magnitude. These reference genes were expressed at comparable levels or lower than the HKGs examined. For accurate determination of inter-assay variations and primer efficiencies, flanking regions of the genes (~300 bp) were amplified by PCR, sub-cloned and the sequences verified. These recombinant plasmids were then linearized and served as templates to construct standard curves. All the qPCR assays showed high efficiency of amplification (>90%) and low intra- and inter-assay variations (Additional File 1). All RNA samples showed RQI values of greater than 9, indicative of high quality and integrity (data not shown).

Stabilities of candidate reference genes and common housekeeping genes

Using both geNorm and NormFinder, we analyzed the expression stabilities of the twenty candidate reference genes and the two commonly used HKGs across all six differentiation conditions. Both statistical approaches recommended the same three ribosomal protein genes RPL19, RPL29 and RPL3 as the overall best reference genes (Figure 3). Pairwise variation analysis by geNorm showed that the combination of RPL19 and RPL29 is sufficiently stable (V2/3 = 0.107, less than the recommended cut-off of 0.15), thus excluding the need to incorporate a third reference gene RPL3 for normalization of target gene expression. Notably, neither GAPDH nor ACTB were recommended.

Further analysis of candidate gene stabilities in each treatment group (Additional File 2A) or at specific time point (Additional File 2B) revealed that the stability rankings of candidate genes do vary among different subgroups. However, with the exception of RPL29 in KCl treated samples, the two genes RPL19 and RPL29 were consistently ranked among top 5 in all subgroups. In contrast, the stability rankings of GAPDH and ACTB varied considerably among different subgroups and they were ranked among the least stable ones within the group of 22 genes in several subgroups. The data indicated that the two novel candidate genes RPL19 and RPL29 have higher expression stabilities than both GAPDH and ACTB, and may serve as better normalizers for gene expression in neuronal differentiation of PC12 cells.

Comparison of the normalization factors generated by different reference gene(s)

To account for possible variations introduced during sample preparation and measurements, raw expression profiles of target genes were scaled by a normalization factor (NF) calculated based on independent measurement of one or more internal reference genes. The variation between NFs generated by different reference genes is thus directly reflective of the variation in the final target gene expression values normalized by different reference genes. We noticed that although RPL19 and RPL29 were ranked as the overall best pair of reference genes, they were not necessarily the best pair for each treatment subgroup. To test the robustness of these two genes across different treatments, we compared the normalization factors calculated based on RPL19 and RPL29 (NF_RPL19/RPL29) to that of the most stable pair of reference genes (NF_top2) in each treatment subgroup. Similarly, we examined the differences between NF_top2 and NFs calculated based on the commonly used HKGs, ACTB (NF_ACTB) or GAPDH (NF_GAPDH). The deviations of each NF from NF_top2 are represented in Figure 4 (see Additional File 3 for details of calculation). The NF_RPL19/RPL29 values were found to least deviate from NF_top2 in NGF, GDNF and KCl subgroups, and had zero deviation in Fsk and Y27632 subgroups as RPL19 and RPL29 were ranked top 2. In contrast, NF_ACTB and NF_GAPDH differed substantially from NF_top2 in many instances, reflective of their varying stabilities across different treatments.

Effect of different reference genes on the interpretation of target gene regulation

Next, the possibility that using scaling factors of NF_RPL19/29, NF_ACTB or NF_GAPDH may substantially alter the interpretation of target gene expression regulation in NGF induced neuronal differentiation was investigated. The relative fold changes of EGR1, ITGA1 and CRYAB expressions normalized by the three NFs were compared to the values normalized by NF of the top 2 genes (NF_RPL29/RPL10A). No statistically significant differences were observed among NF_RPL29/RPL10A, NF_RPL19/RPL29, and NF_ACTB normalized values; whereas NF_GAPDH normalized fold changes were significantly different (Figure 5A-C). In the case of EGR1 and ITGA1, the use of GAPDH as reference gene resulted in the underestimation of target genes expressions, leading to false negative conclusions when a two-fold cut off was applied (Figure 5A-B). On the other hand, normalization by GAPDH resulted in the significant over-estimation of the down-regulation of CRYAB in NGF treated samples (Figure 5C).

The clearly different expression profiles of EGR1, ITGA1, and CRYAB when normalized to GAPDH raised the possibility that GAPDH expression could be regulated over the course of NGF induced differentiation. Normalization of GAPDH expression by the NF of the top 2 genes (NF_RPL29/RPL10A) and the NF_RPL19/RPL29, revealed that GAPDH expression was indeed significantly elevated (>2.5 fold at 24 h) in NGF-stimulated PC12 cells (Figure 6). A more detailed analysis of the kinetics of GAPDH expression over time revealed that expression of GAPDH indeed increased over a period of 28 h (Additional File 4). As a result, the use of GAPDH as a single, unverified reference gene would invariably lead to erroneous interpretation of target gene regulation.

Similarly, we investigated the effect of different reference gene(s) on normalized target gene expressions in GDNF, Forskolin, KCl and ROCK inhibitor Y27632 treated samples. Similar to the case of NGF treatment, with GDNF stimulated PC12-GFRa1a/RET9 and PC12-GFRa1a/RET51 cells, normalization by GAPDH resulted in the underestimations of GDNF induced upregulation of EGR1 and ITGA1 expressions; and over-estimated CRYAB down-regulation (Additional File 5A, B). Interestingly, normalization by ACTB was found to overestimate the expression of EGR1 and ITGA1 expressions in PC12-GFRa1a/RET51 but not RET9 cells (Additional File 5B, 72 h), highlighting the subtle differences between GFRa1a/RET9 and GFRa1a/RET51 systems.

For Forskolin and ROCK inhibitor Y23672 differentiated samples, normalization by ACTB consistently led to the over-estimations of target gene expressions (Additional File 5C, D). Depending on the time point analyzed, normalization by GAPDH was shown to result in either underestimations or overestimations of target gene expressions (Additional File 5C, D). In KCl treated samples, no statistical significant difference was observed among NF_{top2 (RPL19/REPS1)}, NF_RPL19/RPL29, NF_ACTB and NF_GAPDH normalized target gene expression, which suggested that all four were acceptable reference gene(s) for this particular experimental condition (data not shown).

It is thus evident that the stabilities of the two most commonly used HKGs, GAPDH and ACTB vary across different experimental conditions during neuronal differentiation of PC12 cells. They were acceptable reference genes under some conditions but may significantly under- or over-estimate target gene expression under others. On the contrary, the two novel candidate reference genes RPL19 and RPL29 were stably expressed among all conditions analyzed and allowed accurate normalization of differentially regulated genes during PC12 differentiation. It is worthy to note that at early time points (0.5 h and 6 h), the expressions of EGR1, ITGA1 and CRYAB did not show any significant differences when scaled with either NF_top2, NF_RPL19/RPL29, NF_ACTB or NF_GAPDH (data not shown). This is consistent with the observation that the expression of GAPDH did not change significantly at the early time points (Figure 6. and Additional File 4).

Cell Culture

The rat pheochromocytoma cell line PC12 (catalog # CRL-1721; American Type Culture Collection) cells were grown in DMEM supplemented with 10% heat-inactivated fetal bovine serum (FBS; Hyclone, Logan, UT) and 5% Horse Serum (HS), in a humidified atmosphere with 5% CO₂ at 37°C. Wild type PC12 cells, that do not endogenously express GFRa or RET (data not shown), were co-infected with murine GFRa1a (NM_010279) and RET9 (NM_001080780) or RET51 (NM_009050) in pQCXIN or pQCXIH vector by retro-viral infection (Clontech, Mountain View, CA) and selected with 0.4 mg/ml G418 and 0.1 mg/ml Hygromycin, over a period of 2 months.

Differentiation, sample collection and assessment of neurite outgrowth

Two million wild type or infected PC12 cells were seeded in 75 cm² flask (NUNC, Finland) overnight, in DMEM supplemented with 10% FBS and 5% HS, followed by serum depletion for 12 h. PC12 cells were then treated with 50 ng/ml recombinant human GDNF (Peprotech, NJ), 50 ng NGF (Peprotech, NJ), 10 μM Forskolin (Sigma, St. Louis, MO), 5 mM KCl (Sigma, St. Louis, MO) or 25 μM ROCK inhibitor Y27632 (Calbiochem, USA) in DMEM to induce neurite outgrowth. Total RNA was isolated from control and treated cells at 0.5, 6, 24 or 72 h. For neurite outgrowth assessment, cells bearing at least one neurite with the length equivalent to the cell bodies were scored at 72 h by independent observers. More than 400 cells from three different fields were counted per flask.

Immunocytochemistry

Control and treated PC12 Cells were fixed with 4% paraformaldehyde in 1xPBS for 15 min at 37°C, and subsequently in methanol at -20°C for an additional 15 min. After three washes with 1xPBS, cells were permeabilized and blocked with normal goat serum (1:10; Dako, Glostrup, Denmark) in 0.5% Triton X-100/PBS for 30 min at 37°C. The cells were then incubated with high-molecular-weight neurofilament protein (NF-200) antibody (Sigma, St. Louis, MO) at 1:80 dilution in 0.1% Triton X-100/0.1% BSA/1xPBS for 1 h at 37°C and washed three times in 1xPBS. Subsequently, the cells were incubated with goat anti-rabbit fluorescent secondary antibody (Alexa Fluor 488; Invitrogen, CA) diluted 1:1000 in 0.1% Triton X-100/0.1% BSA/1xPBS for 1 h. The cells were then washed three times in 1xPBS. Image acquisition was performed using the Zeiss inverted Axovert 25 microscope equipped with fluorescence detection (Oberkochen, Germany).

RNA Purification and cDNA Preparation

Total RNA from PC12 cells was prepared using TRIzol reagent (Invitrogen, CA) according to manufacturer's instruction. Total RNA was collected from samples in quadruplicate at each treatment time point and the integrity of the RNA validated by denaturing agarose gel electrophoresis and using the StdSens analysis chip on the Experion Automated Electrophoresis System (BioRad, CA) according to manufacturer's instructions. The Experion Automated Electrophoresis System assigns a RQI to each RNA electropherogram which ranges from 10 (intact RNA) to 1 (completely degraded RNA). RNA concentration was quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE), and the 260/280 and 260/230 ratios were examined for protein and solvent contamination. Five micrograms of total RNA were reverse transcribed in a total volume of 20 μl containing 400 U of ImpromII and 0.5 μg of random hexamer (Promega, Madison, WI) for 60 min at 42°C according to the manufacturer's instructions. The reaction was terminated by heating at 70°C for 5 min, and the cDNA was diluted 1:20 for quantitative real-time PCR.

Microarray

PC12 cells were seeded on 25 cm² flask in complete medium and subsequently incubated for 12 h in serum free DMEM. The cells were then treated with GDNF (50 ng/ml) or NGF (50 ng/ml) for 0.5 h or 72 h in duplicates. Total RNA was isolated, quantified and integrity verified before it was amplified using Ambion Illumina RNA Amplification kit (Ambion, TX, USA). Briefly, total RNA (500 ng) was reverse transcribed by ArrayScript in the presence of T7 Oligo(dT) primer. Second strand of the cDNA was synthesized by DNA polymerase at 16°C for 2 h. The cDNA was purified and in vitro transcribed with T7 RNA polymerase and biotin-NTPs. Biotin-labeled cRNA samples were purified and quantified by ND-1000 spectrophotometer (NanoDrop, Fisher Thermo, DE, USA). Each cRNA (750 ng) was hybridized to RatRef-12 Expression BeadChip (Illumina, San Diego, CA, USA) containing 22,523 probes for a total of 21,910 rat genes selected primarily from the NCBI RefSeq database (Release 16) according to instruction provided by Experienced User Card (11286340 Rev A, Illumina). After hybridization, washing and blocking, the BeadChip was incubated with Streptavidin-Cy3 solution (Amersham Biosciences, Piscataway, NJ, USA). Fluorescent signals were obtained from scans on the high resolution Illumina BeadArray reader, using a two-channel, 0.8 μm resolution confocal laser scanner. The Illumina BeadStudio software (Version 2.0) was used to extract fluorescence intensities and the raw fluorescent data was background subtracted and used for analysis. Background is defined as the average signal intensity estimated from the negative control bead types. Outliers are removed using the median absolute deviation method. Detection p-values produced by the BeadStudio software were corrected for multiple hypothesis testing.

Primer Design and Plasmid Standards

The Genbank accession for each candidate reference gene was retrieved from the Illumina microarray probe set and compared to the NCBI RefSeq database (Release 16; http://www.ncbi.nlm.nih.gov. Transcript splicing sites were retrieved from Ensembl http://www.ensembl.org. Where more than one transcript matched the probe, the sequences were aligned and the primers were designed to amplify the consensus region. Vector NTI Advance 10 was used to design two sets of primers for each target gene. The first set of primers generates an amplicon of ~300 bp and is used as a template for RT-qPCR of the targeted gene. The template was subcloned into pGEMT-easy (Promega) vector as previously described [54]. The second set of primers was used for RT-qPCR and was designed to amplify a ~100 bp region within the ~300 bp template. Both primer sets were exon spanning to avoid amplification from genomic sequences. Where possible, primers for RT-qPCR were designed to target the same exons used in the Illumina Expression BeadChip. All primer sequences were evaluated for possible false priming to known rat sequences using the NCBI BLAST tool http://blast.ncbi.nlm.nih.gov/. All products generated after amplifications were verified by gel-electrophoresis and DNA sequencing.

Quantitative Real-Time PCR

Real-time PCR was performed on Biorad iCycler 4 Real-Time PCR Detection System (Bio-Rad, Hercules, CA) using SYBR Green I. The threshold cycles (Ct) were calculated using the iQ5 Optical system software version 2.0. Real-time PCR was performed after an initial denaturation for 10 min at 95°C, followed by 40 cycles of 30 s denaturation at 95°C, 30 s annealing at 60°C, and 30 s extension at 72°C. Fluorescent detection was performed at the annealing phase. Melt curve analysis was carried out at the end of the cycling to confirm that a single product had been amplified. Primer dimer formation in all the assays showed distinct melt characteristics from the correct amplicons. The reaction was performed in a total volume of 40 μl in 1× XtensaMix-SG (BioWORKS, Singapore), containing 2.5 mM MgCl₂, 10 pmol of each primer, and 0.5 U of KlearTaq Hotstart DNA polymerase (KBioscience, UK). All real-time PCR quantification was performed simultaneously with linearized plasmid standards and a non-template control [54]. As PCR is an exponential process, it can be described by the equation, N_n = N₀(1 + ε)ⁿ, where N_n is the number of target molecules at cycle n, N₀ is the initial number of target molecules, ε is the efficiency of amplification and n is the number of cycles. The efficiency of target amplification of an assay was determined from the slope of a plot of Ct (Threshold cycle) versus -log₁₀ concentration of the initial number of target molecules. High efficiency of amplification has a slope approaching the value of 3.32 cycles for every 10-fold dilution of the target. The gene expression levels were interpolated from standard curves and expressed as absolute copy numbers. All Real-time PCR experiments were compliant with the MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments) guidelines [55].

Statistical data analysis

Gene expression stability analysis using two publicly available software tools, geNorm http://medgen.ugent.be/genorm/ and NormFinder http://www.mdl.dk/, were carried out according to authors' instruction.