Cloning and evaluation of reference genes for quantitative real-time PCR analysis in Amorphophallus

Plant materials and treatments

Two types of cultivated species (Amorphophallus albus and Amorphophallus konjac) with a high konjac glucomannan (KGM) content and good quality from the Xiema resource nursery (Xiema, Beibei, Chongqing, PR China) were used as experimental materials. The corms were planted in pots containing soil and peat (2:1, V/V), cultured in a greenhouse and irrigated once every 7 days. After the blades were fully extended, the plants were treated with a simulated high temperature (40 °C) or waterlogging (2.5 cm above the soil surface), and the leaves were harvested after 0 h, 2 h, 4 h, 8 h, and 1 d of exposure (Gu et al., 2011). Different tissue samples (leaves, roots, and corms) were harvested from the same plant after being excavated from the pot and rinsed with distilled water. The materials for verifying the expression profile of SHSP were treated at 42 °C and harvested at 0 h, 1 h, 2 h, 4 h, and 8 h. All samples were immediately frozen in liquid nitrogen after harvest and stored at −80 °C prior to use.

RNA isolation and first-strand cDNA synthesis

Total RNA was extracted from Amorphophallus leaves using the ZH120 Quick RNA Islation Kit (Waryong, Beijing, China). RNA was treated with RNase-free DNase I (Waryong, Beijing, China). Total RNA from roots and ground corms was extracted by the pine method (Chang, Puryear & Cairney, 1993; Gille et al., 2011). The concentration and contaminants were measured using a NanoDropTM 2000 spectrophotometer (Thermo Scientific, USA), and only samples with an OD260/280 between 1.8 and 2.0 and OD260/230 >2.0 were used for cDNA synthesis. The quality and integrity of total RNA were verified by 1% agarose gel electrophoresis. First-strand cDNA synthesis was conducted by using a PrimeScript ™ RT reagent Kit RR047A with gDNA Eraser (TaKaRa, Japan) and generated by using random hexamers. One microgram of total RNA from each sample was used in reverse transcription reactions according to the manufacturer’s instructions. The synthesized cDNAs were diluted to 25-fold for RT-qRCR analyses.

Gene cloning

Degenerate primers of 8 candidate reference genes were designed based on conserved regions of homologous sequences of different monocots, such as Anthurium andraeanum, Phoenix dactylifera, Elaeis guineensis, and Musa acuminata. The accession numbers from GenBank of 4 monocots are shown in Table 1. The accession numbers of the primer sequences are listed in Table 2. The amplified PCR products were A-T cloned into the pEASY-T1 simple vector, and then, the products were transformed into T1 competent cells. The pEASY-T1 simple vector and T1 competent cells were used with a pEASY-T1 simple cloning kit (Transgen, Beijing, China). Positive colonies were identified using colony PCR. The M13 forward and reverse primers were used to recover polymorphisms of the amplification results. Bacterial liquid cultures of positive clone products were sequenced by the Beijing Genomics Institute (BGI, China), and the sequences were identified as the candidate homologous reference genes using NCBI BLAST (Table 3). All of the sequences of the 8 candidate reference genes as well as the published RP sequences (254998327) are shown in Sequence S1.

RT-qPCR assay

RT-qPCR reactions were executed in 96-well plates using the Bio-Rad CFX96 Real-Time PCR system (Bio-Rad, USA) and SsoFast EvaGreen Supermix (Bio-Rad, USA). The total volume comprised 4 ng of cDNA template, 0.2 μM reverse primer, 0.2 μM forward primer, 5 μL of SYBR Green mix, and ddH₂O to 10 µL. RT-qPCR reactions were conducted using the following parameters: 95 °C for 30 s followed by 40 cycles of 95 °C for 5 s and 55−65 °C for 5 s. At the end of the process, the specificity of the amplification was tested using melting curve analysis from 65 to 95 °C with a 0.5 °C increase in temperature at each step. Each RT-qPCR reaction set included water as a negative no-template control for each primer pair. To obtain more reliable experimental results, each PCR assay was conducted in triplicate (technical replicates).

Data analysis

The methods commonly used to analyze the stability of reference genes are geNorm (Vandesompele et al., 2002), BestKeeper (Pfaffl et al., 2004), and NormFinder (Andersen, Jensen & Orntoft, 2004), which are based on different computer algorithms to rank candidate reference genes under different experimental conditions (De Spiegelaere et al., 2015). The geNorm, BestKeeper, and NormFinder programs that were used to calculate the stable values were downloaded for experimental analysis (e.g., geNorm version 3.5 (http://medgen.ugent.be/ jvdesomp/genorm/), BestKeeper version 1 (http://www.gene-quantification.de/bestkeeper.html), NormFinder version 0.953 (http://www.mdl.dk/publicationsnormfinder.htm)).

The expression levels of the candidate reference genes were measured according to their quantification cycle (Cq) values. For geNorm and NormFinder, the raw Cq values were converted to the required data using the formula: 2^−△Cq (△Cq = each corresponding Cq value-minimum Cq value; Liu et al., 2012). The geNorm algorithm was used to further determine the expression stability of the candidate reference genes by calculating the average expression stability value (M) and pairwise variation (V) between all pairs of genes. The ranking of candidate reference genes is based on those with lower M values. In addition, the BestKeeper algorithm used the untransformed Cq values to analyze stability.

Normalization of SHSP gene

The gene expression levels of SHSP were quantified at different time points under heat treatment using the best four reference genes (EF1-a, EIF4A, H3, and, UBQ) separately and a combination of multiple reference genes (EF1-a+H3, EF1-a+UBQ, and EF1-a+H3 +UBQ). Two traditional reference genes (ACTB and GADPH) were also selected for normalization to control for different normalization factors. The software qBasePlus version 3.0 (http://www.biogazelle.com/) was downloaded for these calculations (Hellemans et al., 2007).

Verification of the primer specificity and PCR amplification efficiency

Specific primers for the reference genes and SHSP were designed using Primer3Plus (Untergasser et al., 2007; http://primer3plus.com/) based on the results of sequencing and the published ribosomal protein (RP) sequences of Amorphophallus. The product sizes of the 9 candidate reference genes and SHSP amplified by specific primers were between 100 and 200 bp in length. Each primer pair was evaluated by the presence of a single peak in melting curve analysis (Fig. S1A) and a single, distinct band on a 3% agarose gel (Fig. S1B). The gene-specific PCR amplification efficiency (E) and correlation coefficient (R²) were calculated using a standard curve in which another replicate was performed using a standard curve generated by 10-fold serial dilutions of gel-extracted PCR products. The amplification efficiency (E) was calculated as follows: E = (10^−1∕slope − 1) × 100% (Pfaffl, 2001). Only primers with an ideal value range (110% ≥ E ≥ 90%) and correlation coefficient (R² ≥ 0.99) were used for subsequent experiments. The primer sequences and amplification characteristics, including the Tm, length, efficiency, and R², are shown in Table 4. The PCR amplification efficiencies (E) of nine candidate reference genes and SHSP were between 92.8% and 109.5%. Standard curve regression equation correlation coefficients (R²) were between 0.990 and 0.999, which demonstrated a strong linear relationship.

Cq values of candidate reference genes

The raw Cq values of 9 candidate reference genes ranged from 18.73 (RP) to 34.78 (ACTB; Table S1), and the mean Cq value was used for further analysis. The mean Cq values of the reference genes were between 22.78 (RP) and 27.96 (TUB). The threshold fluorescence for TUB was slightly higher than that of the other genes, indicating that TUB had a low level of expression. In general, there were large variations in the expression levels of each of the 9 reference genes (Fig. 1).

geNorm analysis

geNorm is a Microsoft-based VBA macro software that was developed by Vandesompele et al. (2002). The basic principle is that the ratio of the expression levels of two ideal reference genes under any experimental conditions or cells should be identical in all samples, keeping the calculated M value of a single gene as low as possible. In addition, an increasing number of studies have found that using two or more reference genes contributes to correct system deviation and obtaining more reliable results (Bustin et al., 2009; Vandesompele et al., 2002). Although geNorm can use standardized factors to pair differences, obtaining a lower threshold is preferable.

We ranked the 9 candidate reference genes in 7 sets according to their expression stability values from low to high (Table 5). The threshold of the M value was 1.5; genes with an M value under 1.5 could be considered to be reference genes. EF1-a and EIF4A were two of the best three reference genes in all sets, with values under the threshold in each case. In A. konjac heat-treated samples, the values of all of the candidate reference genes were under the threshold value, and H3 had the lowest value (M = 0.78). UBQ performed well under heat stress in all heat-treated samples. Under waterlogged conditions, TUB was the highest-ranked gene in the two species besides EF1-a and EIF4A, and 9 and 5 candidate reference genes could be regarded as reference genes based on their M values in A. albus and A. konjac, respectively. In a comparison of the expression of the 9 candidate reference genes across different tissues, only EF1-a and EIF4A were under the threshold value in the two species of Amorphophallus. Although UBQ was the best candidate reference gene for A. albus, it was unstable in A. konjac. The M value of CYP was above 1.5 in A. konjac, but it ranked highly in the two species across different tissues. Therefore, we could regard it as an available additional reference gene across different tissues. “Total” contained all test samples, and EIF4A and EF1-a were the best two reference genes in the total set. The values of UBQ and H3 were also under 1.5. The two most stable genes calculated by geNorm at each step during stepwise exclusion of the least stable reference gene are shown in Fig. 2. Starting from the least stable gene at the left, the genes are ranked according to increasing expression stability, ending with the two most stable genes on the right. In general, EIF4A and EF1-a showed remarkable stability in all sets. At the same time, H3 and UBQ under heat stress, TUB under waterlogging, and CYP across different tissues were also expressed stably. In order to avoid using different genes belonging to the same biological process as reference genes, we excluded EIF4A and repeated the geNorm analysis. The results showed that EF1-a was still one of the stable ones in all sets (Fig. 3).

Furthermore, to determine the optimal number of reference genes required for effective normalization, pairwise variation (Vn/Vn+1), which is required for normalization between sequential normalization factors (NFs), was introduced by geNorm. All of the results of pairwise variation are illustrated in Fig. 4. The cut-off value set by the algorithm is 0.15, below which the inclusion of an additional reference gene is not required. For example, V2/3 of A. konjac in heat stress and V4/5 in A. albus under waterlogging were under the cut-off values, which indicated that 2 and 4 reference genes were required for more reliable normalization in these two conditions, respectively. While the 0.15 threshold was not met in the other samples, the geNorm developers emphasized in the user manual, that the proposed threshold of 0.15 must not be taken as a strict cutoff. The cut-off value was set only to offer guidance for determining the optimal number of reference genes. Therefore, they recommend that using only the 3 best reference genes is, in most cases, a valid normalization strategy and results in a much more accurate and reliable normalization compared to the use of only a single reference gene.

NormFinder analysis

Claus et al. compiled the NormFinder program in 2004, and the principle of the program is similar to that of geNorm. NormFinder generates a stable value for gene expression, then sorts gene expression in ascending order according to the stable value (Andersen, Jensen & Orntoft, 2004). The gene with the lowest stable value is the most stable reference gene. However, this program has the drawback that it identifies only a single best reference gene.

The results calculated by NormFinder were similar to those calculated by geNorm. All of the values, ranked from low to high, are shown in Table 6. The lower the value, the higher the stability. In A. albus under heat stress, EF1-a was the best reference gene, followed by EIF4A. In heat-treated A. konjac, H3 was the best reference gene, followed also by EIF4A. Under waterlogging conditions, the results showed that EF1-a and EIF4A were the best two reference genes in the two species. Additionally, UBQ and RP performed well across different tissues in A. albus based on the results calculated by NormFinder. In the total set, EIF4A and EF1-a showed remarkable expression stability in the two species of Amorphophallus. Overall, EIF4A and EF1-a performed very well in all sets and were identified as the best two reference genes in 5 sets.

BestKeeper analysis

BestKeeper software was written by Pfaffl et al. (2004). The data were entered into the BestKeeper Excel file. Then, BestKeeper calculated the standard deviation (SD) and coefficient of variance (CV) of each gene. The candidate reference gene with the lowest coefficient of variance and standard deviation (CV ± SD) was considered to be the best reference gene. Any reference gene with a SD1 was excluded because gene expression was not consistent in all samples. The advantage of BestKeeper software is that it is not only able to analyze the stability of reference genes, but it can also compare the expression levels of target genes.

The results of BestKeeper analysis were different from those of the other two programs which many researches had drew similar conclusions (Lin et al., 2014; Qi et al., 2016). This may be because the principle of this algorithm differs from that of the others. Lower CV values represent higher stability. EIF4A was the best reference gene under heat stress in the two species of Amorphophallus, and ACTB and EF1-α ranked second in A. albus and A. konjac, respectively. In the waterlogging sets, even though EIF4A performed well in A. albus, it was identified as the worst reference gene in A. konjac. At the same time, EF1-α was one of the least stable reference genes under waterlogging in both species of Amorphophallus. Across different tissues, CYP, EF1-α, and EIF4A were ranked the same and were the best three genes in both species. In total, H3, TUB, and UBQ were ranked above EF1-α and EIF4A, even though EF1-α and EIF4A were ranked as unstable in other sets and were even the worst two genes in the waterlogging treatment of A. konjac. The stability of these two genes could be better than others when SD ≤ 1 was taken into consideration. In fact, this algorithm does not consider internal differences between different plants. The calculated results are only references for selecting reference genes. The ranking is shown in Table 7.

Evaluation by normalizing of SHSP

Small heat shock proteins (SHSP) are important proteins in plants that function as molecular chaperones (Jakob et al., 1993). SHSP genes in plants are involved in resistance to abiotic stress and are rapidly induced under heat stress (Sun, Van Montagu & Verbruggen, 2002). The sequences of a SHSP gene from two species were cloned and are shown in Sequence S1. The relative expression levels of SHSP analyzed using 9 different normalization factors are shown in Fig. 5. The results revealed that the trends of SHSP expression were similar by using different normalization factors. However, differences were observed at 1 h, 2 h, and 8 h and the relative expression using the identified reference genes showed almost identical results. The relative expression patterns of normalized by EIF4A and EF1-a were more similar than UBQ and H3. When ACTB and GADPH were used as reference genes, the expression of SHSP at these time points were significantly different from the results using identified stable genes in two species based on the results of student’s T test. Lower sum of squares of deviations (SS) value indicated the data was closer to the expected. The result showed SS values of 3 normalization combinations were lower than these of single identified stable reference genes.