Selection and validation of reference genes for quantitative gene expression studies in Erythroxylum coca

Plant material

Erythroxylum coca was obtained from the Bonn Botanical Garden. Plants were grown at 22°C under a photoperiod of 12 h light/12 h dark with relative humidities of 65% and 70% for light and dark conditions respectively (and fertilized once a week with Ferty 3 (15-10-15) and Wuxal Top N (Planta, Regenstauf, Germany).

The organs used for RNA extraction and qRT-PCR analysis were obtained from four-month old E. coca plants grown from rooted cuttings. Leaves in four developmental stages, roots, stems and flowers were analysed. The leaf developmental stages were: leaf buds; young expanding leaves in a rolled state (Stage 1); young expanded (unrolled) leaves (Stage 2); and fully mature leaves (Stage 3) (see Figure 1).

Figure 1. Developmental leaf stages of Erythroxylum coca plant.

Leaf Stage I (L1) young rolled leaves, Leaf Stage II (L2) young expanded leaves, Leaf Stage III (L3) fully mature leaves.

RNA extraction and cDNA synthesis

Total RNA was extracted from 100 mg of fresh plant tissue using a total RNA extraction kit (Invitek, Berlin, Germany). Genomic DNA was removed by treatment with RNAse-free DNAse I (Qiagen, Hilden, Germany). RNA quality was assessed on an Agilent Bioanalyzer 2100 using a RNA 6000 Nano Kit (Agilent, Böblingen, Germany). RNA concentration was determined using a NanoDrop 2000 c spectrophotometer (NanoDrop Technologies, Wilmington, USA). cDNA was synthesized using a Super Script III First Strand Kit (Invitrogen, Karlsruhe, Germany) according to the manufacturer’s instructions. In brief, random hexamer primers and deoxyribonucleoside5’-triphosphates (dNTPs) were added to 5 µg total RNA and the mixture was incubated at 65°C for 5 min before brief chilling on ice. The first strand was then reverse transcribed by adding First Strand Buffer, 20 mM dithiothreitol and Super Script III reverse transcriptase to a final volume of 20 µl and incubating the mixture at 42°C for 1h. The resulting cDNA was diluted 1:20 (vol:vol) with deionized water and stored at -20°C.

Reference gene selection

Candidate reference genes were selected from an E. coca 454 sequence library² based on their homology to previously reported reference genes in A. thaliana⁹. Nine candidate reference genes with an E-value higher than 2e^-72 were identified by BlastN comparison as orthologues to Arabidopsis genes: Expressed protein (Ec6409), Expressed protein (Ec10131), Clathrin adaptor complex subunit (Ec11142), Actin (ACT), Adenine phosphoribosyl transferase 2 (APT2), Elongation factor 1 alpha (EF1α), Protein tyrosine phosphatase 1B (TPB1), Peroxin 4 (Pex4) and Pp2aa3-like protein. Primers for qRT-PCR were designed using Primer Express 3.0 (Applied Biosystems) and their sequences are shown in Table 1.

All primer pairs were validated prior to their use in gene expression analysis. PCR reactions were performed with each primer pair and the products were visualised by gel electrophoresis to confirm the presence of a single PCR product of the expected size. The sequence specificity of the PCR products was also verified by sequencing.

Quantitative real-time PCR

All PCR reactions were performed on a Stratagene Mx3000P (La Jolla, USA). Each reaction contained 12.5 µl Brilliant Sybr Green (Agilent/Böblingen, Germany), 0.375 µl Rox, 0.4 µM primers and 1 µl cDNA in a final volume of 20 µl. All samples were run in triplicate. The thermocycling conditions were denaturation at 95°C for 10 min; followed by 40 cycles of denaturation (95°C, 15 s) and annealing/extension (60°C, 1 min). A melting curve analysis protocol was performed after completion of the PCR reaction to confirm the absence of multiple amplicons and/or primer dimers. A no template control (NTC) was included to ensure the absence of contamination. In addition, the presence of genomic DNA contamination was excluded by performing reactions without reverse transcriptase. PCR efficiency was determined using a standard curve based on between five and seven different four-fold dilutions of a cDNA cloned amplicon.

Data analysis

Cycle threshold (Ct) values were exported from the MxPro software (Stratagene) to Microsoft Excel using the qBASE v1.3.5 macro¹⁹. PCR efficiencies and regression coefficients were calculated in qBASE and are reported in Table 1. The expression stability of the nine reference genes in E. coca tissues was evaluated with geNorm v3.5¹⁵, NormFinder¹⁶ and Bestkeeeper v1¹⁷. Relative expression quantities were exported from qBASE and analyzed in Microsoft Excel using the geNorm v3.5 and NormFinder macros. For analysis using the BestKeeper macro, Ct values from the MxPro Software and PCR efficiencies calculated by qBASE were utilized.

Selection and expression profiling of candidate reference genes

A similarity search (BlastN) between previously identified reference genes from A. thaliana⁹ and an E. coca 454 sequence library² was conducted to identify orthologous sequences. Nine E. coca genes with high similarity to A. thaliana were selected and PCR primers targeting these sequences were developed (see Table 1). To confirm the specificity of the primers and identity of the amplicons, RT-PCR was performed on cDNA from four developmental leaf stages, stems, roots and flowers. Primer specificity was investigated by electrophoresis and a single amplicon of the expected size was obtained for each primer pair (Supplementary Figure 1). Sequence analysis of ten cloned amplicons revealed that the amplified fragments were identical to the targeted sequences in the 454 sequence database. All primer pairs achieved amplification in fewer than 35 cycles in all samples, demonstrating that all of the candidate reference genes are expressed at experimentally useful levels. The ΔC_t between samples and no template controls (NTCs) was always greater than five cycles, showing that contamination during the setup of the experiment was negligible²⁰. All RNA samples were tested for contamination with genomic DNA by performing qPCR analysis on negative control reverse transcriptase reactions in which the reverse transcriptase was omitted. No amplification product could be detected in these control reactions.

The gene-specific amplification efficiency was calculated by linear regression analysis of the standard curve and ranged between 79% (Ec10131) and 97% (Actin). The coefficient of correlation (r²) of the linear regression analysis was always greater than 0.986 as shown in Table 1, indicating a linear relationship between C_t values and log-transformed transcript quantities in the range of the standard curve.

To ensure that the primer pairs are specific for the desired sequence in all samples and do not target homologous transcripts in some sample subsets, a melting curve analysis of each sample was performed after PCR amplification (Supplementary Figure 2). A single peak in the melting curve specific for each primer pair was obtained for all samples, and no peak could be observed in the melting curves of the control reactions (NTC and negative control reverse transcription reactions).

Expression stability of candidate reference genes

The expression stability of the candidate genes were evaluated with the geNorm, NormFinder and BestKeeper algorithms (Table 2). C_t values were transformed to relative quantities using qBASE prior to analysis with geNorm and NormFinder, while C_t values and PCR efficiencies were used in BestKeeper. The cDNA samples were considered as either a single, diverse set derived from all organ samples; or as two subsets derived from leaf buds and leaves (leaf buds, Stage 1, Stage 2 and Stage 3 leaves) or mature organs (Stage 3 leaves, flowers, roots and stems).

geNorm calculates the average expression stability value (M) for each candidate gene on the basis of the average pair-wise variation between all genes analyzed. geNorm analysis indicated that Ec10131 and Ec6409 are the most stable candidate reference genes across all of the E. coca tissues tested (Table 2). In the leaf bud/leaf sample subset, Ec10131, TPB1 and Ec6409 were ranked as the three most stable genes (in that order) (Supplementary Table 1), while in the mature organ subset Ec10131 and Ec6409 were again ranked as the most stable. In contrast, Pex4 and APT2 were consistently ranked as the least stable in all sample subsets (Table 2 and Supplementary Table 1 and Supplementary Table 2). The ‘housekeeping’ genes Actin and EF1α were relatively unstable and were ranked at positions six and seven (respectively) in all sample sets.

The optimal number of reference genes required for accurate normalization in the respective sample sets (all samples, leaf bud/leaf and mature tissues) was determined by calculating the mean variation in each normalisation factor (V) and then observing the effect of iterative addition of the next most stable reference gene (V_n/V_n+1) (as detailed in Vandesompele et al. 2002¹⁵). In each case, the two most stable reference genes were sufficient for accurate normalization, since inclusion of a third gene had little impact on the calculation of the normalization factor (V_n/V_n+1 below 0.15).

BestKeeper ranks gene stability by calculating the correlation coefficient (r) between the expression of each candidate gene and the BestKeeper index (BKI; calculated using all genes across all samples). Across all of the samples tested, BestKeeper indicated that Actin (r = 0.784) and Ec6409 (r = 0.768) were the most stable, while Ec10131 was ranked as the least stable (r = 0.638). In the leaf bud/leaf sample subset, Actin (r = 0.869) and APT2 (r = 0.868) had the highest correlation with the BKI, and Ec10131 again showed the lowest correlation (r = 0.385). In the mature organs sample subset, Pex4 and APT2 were strongly correlated with the BestKeeper index (r = 0.767 and r = 0.724, respectively), whereas Ec10131 showed low correlation (r = 0.309) (Supplementary Table 1 and Supplementary Table 2).

To provide a further ranking of gene stability, the results were also evaluated with NormFinder, in which candidate reference genes are ranked according the variance of their expression relative to the expression variance within a defined group of samples¹⁶. Pp2aa3 was the most stably expressed gene with the lowest expression variance (stability value of 0.291), followed by Ec6409 and Ec11142, when all samples were included in the calculation. When the leaf bud/leaf and mature organ subsets of samples were considered, the rankings varied considerably (Supplementary Table 1 and Supplementary Table 2). Actin, APT2 and Pex4 were always ranked as the seventh, eighth and ninth most stable reference genes (respectively), but there was no consistent order of ranking for the other reference genes. The NormFinder rankings were also distinct from the geNorm rankings, although both algorithms identified Actin, APT2 and Pex4 as having the least stable expression profiles.