Selection of reference genes for quantitative real-time PCR analysis in halophytic plant Rhizophora apiculata

Plant materials

In the present study, we collected three month old R. apiculata seedlings located in the west coast of India with the geographical latitude of 15.5256°N and longitude of 73.8753°E, with the permission from the Principal Chief Conservator of Forest, Goa Forest Department, Goa, India. Mangrove species identification was performed based on morphological characteristics using a comparative guide to the mangroves of Goa (Naskar & Mandal, 1999). All seedlings were acclimatized and maintained in half-strength Hoagland solution at a temperature regime of 24–30 °C, 40–50% relative humidity. Various physiological tissues such as leaves, stems, roots and flower samples were collected. To imitate salt stress conditions, young seedlings of R. apiculata were exposed to Hoagland nutrient solution supplemented with 250 mM sodium chloride (NaCl) continuously and leaf samples were harvested at different time-course such as 0, 6, 12 and 24 h.

RNA isolation and cDNA synthesis

Total RNA was extracted using modified cetyl-trimethyl ammonium bromide (CTAB) protocol with 2% polyvinylpyrrolidone and 10% β-mercaptoethanol (Fu et al., 2004). Freshly collected tissues were immediately pulverized into 2 ml of pre-warmed CTAB buffer and incubated at 60 °C for 30 min. The suspension was gently mixed and centrifuged at 14,000 rpm for 10 min at room temperature with an equal volume of chloroform: isoamyl alcohol (24:1). The aqueous phase was transferred to a new tube and RNA was precipitated with a 1/3^rd volume of 8M lithium chloride (LiCl) and incubated at −20 °C for 1 h followed by adding an equal volume of chilled isopropanol (−20 °C). The RNA was precipitated by centrifugation at 14,000 rpm for 10 min at room temperature followed by washing with 70% ethanol. RNA was finally dissolved in 0.1% DEPC treated water and its quantity and quality were confirmed by Nanodrop (Thermo Fisher Scientific, Waltham, MA, USA). Genomic DNA contamination was removed by DNase I enzyme (Thermo Fisher Scientific, Waltham, MA, USA) treatment at 37 °C for 30 min and heat inactivated at 65 °C for 10 min with 50 mM EDTA. The cDNA synthesis was performed in 20 µl reaction volume using the RevertAid Reverse Transcriptase (Thermo Fisher Scientific, Watham, MA, USA), 0.1–5 µg RNA sample and oligo d(T)₁₈ primer, as per manufacturer’s instructions.

Selection of reference genes and primer designing

Nine housekeeping genes such as ACT,α-TUB, β-TUB, GAPDH, UBQ, 18S rRNA, rbcL, Histone H3 and EF1α used in qRT-PCR along with one target gene sodium/proton antiporter (NHX) were selected. There is no genome sequence available publicly for R. apiculata hence, homologous candidate reference gene sequences were retrieved from model plants such as Arabidopsis thaliana and Oryza sativa from Gramene and NCBI databases. Full-length candidate reference gene sequences were used for primer designing using PrimerQuest (Integrated DNA Technologies) with given parameters: melting temperature (T_m) of 55–65 °C, primer length of 17–25 bp, and amplicon length of 100–500 bp (Table 1). The amplicon was sequenced and annotated based on the sequence similarity-based search tool. Further, all the confirmed sequences were submitted to GenBank for accession numbers. After primer specificity analysis α-Tub and Histone H3 were removed from further analysis. The primer sequences, accession numbers, and their efficiency were given in Table 1.

For qRT-PCR, primer specificity was determined using melting curve analysis and the PCR products were checked on 2% agarose gel. The primer efficiency of all candidate reference genes was calculated based on the standard curve generated from a 10-fold serial dilution of cDNA (10⁰, 10⁻¹, 10⁻², and 10⁻³) and regression coefficient (R²) values. Primer efficiency was calculated using the given formula [E = (10^{(−1∕slope)} − 1) × 100], where E = 2 and corresponds to 100% efficiency; high/acceptable amplification efficiency equals 90–110% (Sinha et al., 2015).

Quantitative RT-PCR analysis

The quantitative RT-PCR analysis was carried out using SYBR green master mix (2X Brilliant III SYBR^® Green QPCR; Agilent Technologies, Santa Clara, CA, USA), on AriaMx Agilent system (AriaMx; Agilent Technologies, Santa Clara, CA, USA) with the following reaction conditions: initial denaturation at 95 °C for 3 min, 40 cycles of 95 °C for 30 s, 55 °C for 30 s and 72 °C for 45 s extension, and a melt-curve program (65–95 °C with a temperature increase of 0.5 °C after every 5 s). The melting curve was generated to determine the amplicon specificity. The qRT-PCR experiments were performed using three biological and two technical replicates. A reaction with no template control and a reverse transcription negative control were performed to check the potential reagents and genomic DNA contamination.

Analysis of gene expression stability

The candidate reference gene ranking was analyzed using five different algorithms such as geNorm, NormFinder, Bestkeeper, Δct, and comprehensive ranking analysis by RefFinder.

geNorm analysis

The geNorm determines the most stable reference genes based on the gene expression stability value (M) for a reference gene. It also calculates the minimum number of candidate reference genes required for normalization of target genes. It requires calculated Cq values into relative quantities using the given formula: Q = E^(minCq−Cq), where Q represents sample quantity relative to sample with the highest expression, E is amplification efficiency and min Cq is the lowest Cq values. The stability value (M) is defined as an average pairwise variation (V) of the gene compared with all other tested reference genes and the cut-off is 1.5 (Vandesompele et al., 2002). If M value is lower than 1.5, it represents stable candidate reference gene and higher values reflect least stable.

NormFinder

NormFinder calculates expression stability values for candidate reference genes and evaluates the most stable reference gene pairs. It also calculates intra and intergroup variation using a direct comparison between genes. It uses same input calculation files which are required for geNorm with a little variation such as the first row represents a sample, the first column represents genes and the last row represents a group of samples. NormFinder is available with Excel spreadsheet add-in (https://moma.dk/normfinder-software). It ranks candidate reference genes based on expression stability value. Lowest M value represents most stable reference gene and higher the value, least stable are the genes (Andersen, Jensen & Orntoft, 2004).

BestKeeper analysis

BestKeeper determines the best reference gene based on the normalization factor (also called Bestkeeper index) and pairwise correlation analysis. It requires raw Cq values as an input data to select most stable and least stable candidate reference genes. It is available in an MS Excel spreadsheet file (http://www.gene-quantification.de/bestkeeper.html) and in RefFinder (http://150.216.56.64/referencegene.php?type=reference) as well. It evaluates the candidate reference gene stability by comparing the standard deviation of each gene and averages of these values. It also calculates the coefficient of variance, Pearson correlation coefficient (r) values, geometric mean (GM) and arithmetic mean (AM).

ΔCt method

This tool is available in an MS Excel spreadsheet as well as in RefFinder (http://150.216.56.64/referencegene.php?type=reference). It calculates the stable candidate reference gene based on standard deviation and pairwise comparison with other genes. ΔCt requires raw Cq values as an input data. It considers a pair of gene for calculations and compares ΔCt values among genes (Silver et al., 2006).

RefFinder analysis

RefFinder is a web-based comprehensive tool developed for evaluating and screening reference genes from extensive experimental datasets. RefFinder was used to generate comprehensive stability rankings (Xie et al., 2012). Comprehensive ranking of seven candidate reference genes was analyzed using RefFinder.

Validation of candidate reference genes

The reliability of highly stable candidate reference genes identified in the current study was validated using sodium/proton antiporter (NHX) as a salt stress target gene. The differential gene expression profiles of NHX under salt stress at 0, 6, 12 and 24 h were normalized using EF1α, ACT, 18S and UBQ along with the combination of EF1α + ACT genes. The input values for EF1α + ACT were calculated using the geometric mean formula given below to normalize gene of interest Geometric $Mean=\sqrt[n]{\times 1,}\times 2,\times 3\dots .\times n$, where the n = number of times (Vandesompele et al., 2002). The average Cq values from three biological replicates were used for relative expression analysis and the relative gene expression level calculated using the 2^−ΔΔCT method (Livak & Schmittgen, 2001; Manuka, Saddhe & Kumar, 2018). Statistical analysis was performed using SPSS 15.0 for Windows evaluation version to verify the significant difference between relative gene expressions. One-way Analysis of variance (ANOVA) with Tukey’s Honest Significant Difference (HSD) test was performed for comparison between reference genes and target genes. A p-value <0.05 was considered statistically significant.

Minimum Information for publication of qRT-PCR experiments guidelines (MIQE)

All the qRT-PCR experiments and data analysis in the present study were performed in accordance with the MIQE guidelines (Bustin, 2002).

Expression profiling of selected reference genes

In order to select stable reference genes, transcript levels in tissues such as leaf, root, stem, and flower as well as salt stress samples were quantified based on their cDNA concentration. The primer specificity was determined by PCR products wherein single, expected amplicon size was obtained (Fig. 1). The qRT-PCR melting curve for template test and negative control (NTC) without template were analyzed for primer-dimer and reagents contamination (Fig. S1). Further, NTC samples were confirmed by running 2% agarose gel electrophoresis. The amplified PCR products were sequenced and submitted to GenBank for accession numbers. All the sequenced PCR products were identified and annotated based on BLAST search. The primer efficiency (%) ranged from 103.70%, (R² = 0.997) for NHX to 90.54% (R² = 0.98) for UBQ including 18S (102.30, R² = 1), ACT (92.03%, R² = 0.994), EF1α (98.60%, R² = 0.99), β-TUB (94.08%, R² = 0.996), GAPDH (96.78%, R² = 0.992) and RbcL (97.69%, R² = 0.996) (Fig. S2; Table 1).The mean cycle threshold (Cq) values of the seven selected reference genes for different tissue samples ranged from 14.16 for18S to 21.77 for GAPDH (Fig. 2A). Similarly, for the salt stress samples, the mean Cq values ranged from 13.96 for 18S to 24.23 for UBQ (Fig. 2B). Mean Cq values gave insight into approximate gene expression data. Negative control showed higher Cq values indicating no product amplification which was further checked on a 2% agarose gel. Moreover, negative control without reverse transcriptase did not show any product amplification, thus indicating no gDNA contamination.

geNorm analysis

For physiological tissues, seven candidate reference genes showed average expression stability value (M) less than 1.5. ACT (M = 0.721) was most stable reference gene followed by EF1α (M = 0.761), and β-TUB (M = 0.763) (Table 2A; Fig. S3A). GAPDH was the least stable candidate reference gene with M value 1.599. The geNorm also determines an optimum number of candidate reference gene for normalization based on the calculation of pairwise variation (V_n∕V_n+1) between sequential normalization factor (NF_n and NF_n+1). To select the best pair for normalization, the threshold value is 0.15. If pairwise variation value is lower than 0.15, there is no need to add more candidate reference gene. Moreover, the best pairwise variation value 0.382 was observed for a combination of ACT and EF1α and comprehensively recommended for normalization (Table 2B; Fig. S3A). Based on the observation, there were no effects on an addition of the third gene in the combination of ACT and EF1α which showed pairwise variation value below 0.15 (Fig. 3A).

Under salinity stress, ACT was most stable candidate reference gene with M value 0.927, followed by 18S and EF1α showing same stability M value 1.02 (Table 2A). Moreover, rbcL and UBQ performed least stable candidate reference gene with M value 1.357 and 1.399 respectively. In salt stress, ACT + EF1α were the most suitable combination for normalization of the gene of interest with pairwise variation value of 0.462 (Table 2B; Fig. S3B). According to pairwise variation analysis, if the third gene was added in the ACT + EF1α, it showed higher pairwise variation value of 0.215 (Fig. 3B).

NormFinder

In R. apiculata physiological tissue samples EF1α was most stable with stability value of 0.085. β-TUB (0.135) was the second most stable candidate reference gene followed by ACT (0.164) (Table 3A). EF1α and β-TUB (0.070) showed the most stable combination for the pair of candidate genes for normalization (Table 3B). Overall, GAPDH, UBQ, and RbcL were least stable reference genes. In salt stress, ACT was most stable reference gene with a stability value of 0.196. EF1α and 18S were second and third most stable candidate reference genes with stability value 0.257 and 0.273 respectively. ACT and EF1α showed the best pair of reference genes with stability value 0.183 (Table 3B). Under salt stress, geNorm and Normfinder showed almost similar results for a selection of candidate reference gene.

BestKeeper

In the BestKeeper analysis, standard deviation (SD) and coefficient of correlation (r) value were the criteria used for comparison. Highest r value represents the most stable candidate reference genes and lower r value represents the least stable genes. Here, we considered r value for evaluation, showing EF1α as the most stable reference gene followed by ACT with r value 0.987 and 0.966 respectively. GAPDH was ranked as the least stable candidate reference gene with lower r value (Table 4). The result is consistent with geNorm and NormFinder analysis. In salt stress, ACT (r = 0.638) showed most stability followed by β-TUB (r = 0.625) and EF1α (r = 0.523) (Table 4). Under salt stress, similar results were observed with little variation in BestKeeper. BestKeeper determined β-TUB second most stable candidate reference gene in salt stress.

ΔCt analysis

According to ΔCt analysis, EF1α was the most stable candidate reference gene followed by ACT and β-TUB in physiological tissue (Table 4). 18S was ranked as an average or moderately stable reference gene. The results were consistent with earlier analysis. GAPDH, RbcL, and UBQ were the least stable. Under salt stress, ACT was most stable candidate reference gene followed by EF1α and18S (Table 4).

Comprehensive ranking of candidate reference genes

Based on the Geomean value, a comprehensive ranking of all candidate reference genes showed EF1α (1.32) was the most stable followed by ACT (2.34) and β-TUB (2.91) (Fig. 4A). Moreover, rbcL and GAPDH performed as a least stable candidate reference genes. In salt stress, ACT was the most stable with Geomean value 1. Moreover, 18S was second most stable candidate reference gene with Geomean value 2.29 (Fig. 4B). Here, UBQ and rbcL performed as least stable candidate reference genes.

Validation of stable candidate reference genes under salt stress

To validate the efficacy of candidate reference genes, ACT, EF1α, 18S and UBQ were used to normalize the expression levels of NHX in salt stress at four different time course. Set of the most stable candidate reference genes such as ACT, EF1α, 18S and the least stable candidate reference gene UBQ were used as internal controls. While using EF1α, ACT, and 18S alone for normalization, NHX showed significant upregulation expression pattern in salt stress at 12 h. However, with UBQ as an internal control, NHX expression was upregulated in salt stress after 6 h of salt stress (Fig. 5).