Reference genes for qRT-PCR normalisation in different tissues, developmental stages, and stress conditions of Hypericum perforatum

Plant materials

H. perforatum seeds (2n = 2× = 16) were germinated on a seedling bed in the glasshouse (25 ± 2 °C, natural lighting, 60–80% humidity). Whole plant tissues were collected at the one-month-old (1M), two-month-old (2M), three-month-old (3M), and six-month-old (6M) stages. Samples of different tissues (leaf, flower, stem, and root) were taken from two-year-old plants (2n = 2× = 16). 3M seedlings were subjected to abiotic stresses and hormonal treatments, including 10 μM salicylic acid (SA), 200 μM methyl jasmonate (MeJA), 100 μM abscisic acid (ABA), one mM AgNO₃ (Ag), 200 μM CuSO₄ (Cu), 100 mM NaCl (Na), low temperature (4 °C) and wounding (W). The stress-treated samples were each collected 6 h after the corresponding treatments, and the control groups were collected following non-treatment. All samples were collected in three replicates, frozen in liquid nitrogen immediately, and then stored at −80 °C.

Total RNA isolation and cDNA synthesis

Total RNA was extracted using the Polysaccharide and Polyphenols Plant Quick RNA Isolation Kit (centrifugal column type; Waryong, Beijing, China). The genomic DNA was digested with RNase-free DNase I (TaKaRa, Kusatsu, Japan). The total RNA was quantified using the absorbance at A₂₆₀/A₂₈₀ and A₂₆₀/A₂₃₀ nm measured with a NanoDrop 2000c spectrophotometer (Thermo Scientific, Waltham, MA, USA). A 1% (p/v) agarose gel was run to visualize the integrity of the RNA. Only RNA samples with an A₂₆₀/A₂₈₀ wavelength ratio between 1.9 and 2.1 and an A₂₆₀/A₂₃₀ ratio close to 2.0 were used for cDNA synthesis. Then, 1.0 μg DNA-free total RNA was used to synthesize first-strand cDNAs with a PrimeScript RT Reagent Kit (TaKaRa, Kusatsu, China) in a 20 μL volume. All cDNA samples were diluted (1:40) with DNase/RNase-free deionized water for qRT-PCR.

Selection of reference genes and primer design

We performed genomic sequencing of H. perforatum using Illumina paired-end, 10× Genomics linked reads and PacBio SMART (GenBank accession numbers: MK054303, MK106356–MK106365). The transcriptome sequencing of H. perforatum for the roots, stems, leaves, and flowers assisted annotation. The 14 candidate genes, including nine traditional HKGs (ACT2, ACT3, ACT7, CYP1, EF1-α, GAPDH, TUB-α, TUB-β, and UBC2) and five potential reference genes (PKS1, GSA, RPL13, SAND, and PP2A) which have been used as candidate genes in other studies (He et al., 2016; Li et al., 2017; Velada et al., 2014), were selected for the assessment of the most stably expressed reference genes (Table S1). To ensure the accuracy of the reference gene predictions, we first screened the candidate genes according to the genome annotation of each, which was assigned based on the best match of the alignments using Blastp to SwissProt, KEGG, NR, and TrEMBL databases. Then the coding sequences of the 14 selected genes were used as queries for BLAST orderly through the TAIR database (http://www.arabidopsis.org/) to further ensuring accuracy. The sequences with the highest homology with Arabidopsis are shown in Table 1. The primers of all the genes were designed using GenScript (https://www.genscript.com) with a melting temperature between 59 and 61 °C, a primer length of 20–25 bp, and an amplicon length of 70–180 bp. The descriptions of the candidate reference genes, primer sequences, and qRT-PCR amplification efficiencies are presented in Table 1.

qRT-PCR conditions and analysis

PCR reactions were performed on the Roche LightCycler 96 system using SYBR^® Master Mix. Reactions were performed in triplicate in 20 μL volumes containing five μL 30-fold diluted synthesized cDNA, 10 μL SYBR^® Master Mix, 0.4 μL 10 mM forward primer, 0.4 μL 10 mM reverse primer and 4.2 μL DNase/RNase-free deionized water. The cycling conditions were 95 °C for 30 s, 45 cycles of 95 °C for 5 s, and 60 °C for 30 s, and a final melting curve analysis. Each reaction was set up with a negative control group and triple technical replicates were performed in each of three biological samples. To calculate the amplification efficiency from 10-fold continuous dilutions of the cDNA (10⁰, 10⁻¹, 10⁻², and 10⁻³) for each gene, standard curves were constructed to obtain the correlation coefficients (R²) and slope values. Using these standard curves, the corresponding PCR amplification efficiencies (E) were calculated (E = (10^−1/slope −1) × 100) (Aleksandar et al., 2004).

Assessment of expression stability

The expression stability of the genes was analyzed using three different Visual Basic applets, GeNorm (Schlotter et al., 2009), NormFinder (Obrero et al., 2011; Zhang et al., 2016), and BestKeeper (Rhinn et al., 2008; Rotenberg et al., 2006). GeNorm derives a stability measure (M-value) via the stepwise exclusion of the least stable gene, and creates a stability ranking (Hu et al., 2009). Genes with M < 1.5 are generally considered stable reference genes (Tong et al., 2009). This measure is based on the principle that the expression ratio of two ideal control genes should be identical in all samples; thus, genes with the lowest M-value are the most stably expressed (Dekkers et al., 2012). The GeNorm pairwise variation (V) values between ranked genes (V_n/V_n+1) were determined for the optimal number of reference genes to be used in the quantitative analysis experiments. A cut-off of 0.15 (V_n value) is usually applied (Vandesompele et al., 2002). NormFinder uses an ANOVA-based model to estimate intra- and intergroup variation within tissues or treatments (“groups” in NormFinder terminology) to assess the expression stability (Marchal et al., 2013). For GeNorm and NormFinder, the raw Ct values need to be converted to relative quantities (Q) using the formula Q = 2^−ΔCt, in which ΔCt = each average Ct value−minimum Ct value (Yang et al., 2015). In BestKeeper, the coefficient of variance (CV) and the standard deviation (SD) were calculated using the Ct values, with lower CV and SD values indicating higher stability (Migocka & Papierniak, 2011). The mean standard error and level of statistical significance were calculated using GraphPad Prism 6.0, and the level of statistical significance was assessed using ^∗P < 0.05, ^∗∗P < 0.01, and ^∗∗∗P < 0.001.

Validation of reference gene stability

The relative expression of the HpHYP1 gene in different tissues was measured and standardized by using the most stable and unstable candidate genes as internal references to verify the reliability of the selected genes according to the 2^−ΔΔCt method (Kumar et al., 2011). Three technical replicates were performed for each biological sample.

Primer specificity and expression level analysis of candidate reference genes

The gene names and abbreviations, accession numbers, primer sequences, amplification efficiencies, amplicon sizes, Tm values, and molecular functions are listed in Table 1 and Table S1. The amplification efficiencies listed in Table S2 ranged from 92.5% (RPL13) to 109.5% (TUB-α), the Tm values varied from 81.2 °C (PP2A) to 87.8 °C (CYP1), and the amplicon sizes were between 76 (GAPDH and PKS1) and 170 bp (TUB-β). Furthermore, primer specificities were determined using melting curves (Fig. S1) and triple technical replicates were performed for each of the three biological samples. A single band indicated the correct size of each pair. The average raw CT values of different genes ranged from 18.15 to 32.32 (Table S3). Data were analyzed within experiments and divided into four groups: tissues from two-year-old plants (TS: R, S, L, and F), seedling developmental stages (SG: 1M, 2M, 3M, and 6M), 3M seedlings exposed to abiotic stresses (ST: SA, MeJA, ABA, Cu, Ag, Na, 4 °C and W) and a combination of all experimental conditions (TT). The expression levels of the 14 selected genes in the four groups are shown in Fig. 1. Among them, UBC2 (TS) had the highest Ct value, but GAPDH (SG) had the lowest, indicating their levels of expression.

Stability of candidate reference genes

GeNorm was used to calculate the normalization factor from the geometric mean of the genes to identify the most stably expressed gene. The M-value is defined as the average pairwise variation of a particular gene with all other potential reference genes. The expression stability ranking of the 14 reference genes in the TT samples was arranged as follows: TUB-β > ACT2 > CYP1 > EF1-α > TUB-α > SAND > UBC2 > ACT7 > PP2A > GSA > PKS1 > RPL13 > GAPDH > ACT3. Except for in the TT group, ACT2 and ACT3 were the genes with the lowest and highest M-values, respectively, (Fig. 2; Table S4). The V_2/3 values for almost all of the experimental sets in H. perforatum were lower than the cut-off threshold of 0.15 (Fig. 3), which indicated that the combination of two reference genes (ACT2 and TUB-β, or ACT2 and EF1-a) can accurately standardize these samples. For the TT group, V_2/3 was 0.2039, which indicated that the top three reference genes (ACT2, TUB-β, and CYP1) were required for accurate normalization.

NormFinder with lower values indicating higher stability is another Excel application that ranks the candidate genes based on their minimal combined inter- and intragroup variation of expression. The results (Table 2) indicated that ACT2 was the most stable gene in the SG and ST groups. EF1-α and TUB-β ranked first in the TS and TT groups respectively. ACT3 was at the end of the rankings in all groups calculated using NormFinder. The stability order of the candidate genes calculated using NormFinder was nearly the same as that calculated using GeNorm, although there was a small difference. For example, ACT2 was identified as the most stable reference gene in the TS group with the GeNorm analysis, while its stability rankings placed it third within the NormFinder analysis.

BestKeeper is another Excel-based program and is able to compare the expression levels of up to 10 target genes, each in up to 100 biological samples. The CV and SD values of the 14 candidate genes computed using BestKeeper are represented in Table 3. The smaller the SD and the CV values are, the better the stability of the internal reference genes. The expression of the reference gene was unstable when SD value is greater than one. The results revealed that ACT2 was the most stable genes for the ST and TT samples, while TUB-β ranked first in the TS and SG samples.

Reference gene validation

To test the reliability of the results, the relative expression patterns of the target gene HpHYP1, which belongs to the PR-10 family and is associated with stress control, were evaluated using different internal control genes in the roots, stems, leaves, and flowers. Based on the results of the present study, the top two stable genes (ACT2 and TUB-β) and the most unstable gene (ACT3) were chosen as internal controls. As shown in Fig. 4, when normalized using ACT2 as the reference genes, the transcript abundance of HpHYP1 was upregulated compared with the results in the root samples. When ACT2 and TUB-β (as identified using GeNorm) were used as the internal references, the expression patterns were similar to those obtain with ACT2. When normalization was based on TUB-β alone, the expression level of HpHYP1 was still up-regulated. The only difference was that the expression level in stems was lower than that in leaves, although the difference was not significant. In general, when stable reference genes, including ACT2, TUB-β, and their combination, were used as internal parameters, HpHYP1 had the highest expression level in flowers. When normalized using the less stable gene ACT3, the target genes in all tissues were expressed at lower levels but showed significantly up-regulation (P < 0.001). Thus, the selection of inappropriate reference genes can lead to over- or under-estimation of the relative transcript level, which might lead to a biased result.