Identification of reliable reference genes for quantitative real‐time PCR analysis of the Rhus chinensis Mill. leaf response to temperature changes

Rhus chinensis Mill. (RCM) is the host plant of Galla chinensis, which is valued in traditional medicine. Environmental temperature directly determines the probability of gallnut formation and RCM growth. At present, there is no experiment to systematically analyse the stability of internal reference gene (RG) expression in RCM. In this experiment, leaves that did not form gallnuts were used as the control group, while leaves that formed gallnuts were used as the experimental group. First, we conducted transcriptome experiments on RCM leaves to obtain 45 103 differential genes and functional enrichment annotations between the two groups. On this basis, this experiment established a transcriptional gene change model of leaves in the process of gallnut formation after being bitten by aphids, and RCM reference candidate genes were screened from RNA sequencing (RNA‐seq) data. This study is based on RCM transcriptome data and evaluates the stability of 11 potential reference genes under cold stress (4 °C) and heat stress (34 °C), using three statistical algorithms (geNorm, NormFinder, and BestKeeper). The results show that GAPDH1 + PP2A2/UBQ are stable reference genes under heat stress, while GAPDH1 + ACT are the most stable under cold stress. This study is the first to screen candidate reference genes in RCM and could help guide future molecular studies in this genus.


Rhus chinensis
Mill. (RCM) is the host plant of Galla chinensis, which is valued in traditional medicine. Environmental temperature directly determines the probability of gallnut formation and RCM growth. At present, there is no experiment to systematically analyse the stability of internal reference gene (RG) expression in RCM. In this experiment, leaves that did not form gallnuts were used as the control group, while leaves that formed gallnuts were used as the experimental group. First, we conducted transcriptome experiments on RCM leaves to obtain 45 103 differential genes and functional enrichment annotations between the two groups. On this basis, this experiment established a transcriptional gene change model of leaves in the process of gallnut formation after being bitten by aphids, and RCM reference candidate genes were screened from RNA sequencing (RNA-seq) data. This study is based on RCM transcriptome data and evaluates the stability of 11 potential reference genes under cold stress (4°C) and heat stress (34°C), using three statistical algorithms (geNorm, NormFinder, and BestKeeper). The results show that GAPDH1 + PP2A2/ UBQ are stable reference genes under heat stress, while GAPDH1 + ACT are the most stable under cold stress. This study is the first to screen candidate reference genes in RCM and could help guide future molecular studies in this genus.
Quantitative real-time PCR (qRT-PCR) is one of the most accurate and commonly used methods to detect mRNA abundance and gene expression and is characterized by high technical sensitivity and precision. The introduction of RT-PCR technology has greatly improved and simplified the quantification of nucleic acids and has become a valuable tool for many different disciplines [1]. Reference gene selection and sample preparation are the key issues that define the reliability of obtained qRT-PCR data [2]. According to existing research, many reference genes (RGs) show unstable expression in different experimental organisms, experimental materials, tissues, and organs, following stress treatment, or at different stages of growth and development [3,4]. Therefore, it is necessary to experimentally test genes that are traditionally thought to be stably expressed under different conditions and experimental treatments.
Rhus chinensis Mill. (RCM), also known as sumac, belongs to the Anacardiaceae family and has a long history of culinary use and as a traditional medicine. Asian practitioners of folk medicine have been using sumac bark for a long time, in addition to its leaves, roots, stems, fruit, and especially the galls on the leaves. Studies have shown that sumac tree compounds have strong antiviral, antibacterial, anticancer, antidiarrhoea, and antioxidant activities [5]. Galla chinensis is a gall formed by the parasitic aphid Schlechtendalia chinensis on RCM leaves, that stimulates the proliferation and expansion of leaf tissue cells. Galla chinensis is considered to have preventive and therapeutic effects against diseases such as diarrhoea, dysentery, rectal and bowel cancer, diabetes, septicaemia, and oral diseases. Galls on RCM leaves are rich in gallotannin (50-70%), a type of hydrolysable tannin. The global annual demand for tannins in different industrial areas, including leather processing, textile printing, and mineral separation, exceeds 15 000 tons. The total annual value of tannin related industries is approximately $30 billion [6]. Gallnut formation involves changes in the expression of many genes, which lead to major changes in hormone signalling pathways and secondary metabolite accumulation that result in the accumulation of tannins [7,8]. Gallnut is an insect gall on RCM leaves, that is mainly parasitized by Melaphis chinensis, and the specific mechanism and influence of temperature on its formation have not been published. However, fruit growers in the artificial planting and harvesting base in Cili County have made long-term observations on the effect of temperature on the formation rate of gallnuts, and this information can be shared with us. We can confirm that temperature has an important effect on gallnut formation, and temperature also affects gene expression and ultimately affects plant-insect interactions and product formation. In addition to the interaction between aphids and RCM, temperature has an important effect on gene expression [9,10]; however, evaluation of the stability of RCM RG expression at different temperatures has not been reported, which hampers research on the changes in gene expression related to gallnut formation and tannin accumulation.
Therefore, in view of the important role of temperature in gallnut formation, the value of gallnut as a traditional medicine and the value of tannins, which are the main components of gallnut, in industry, and other fields, it is extremely important to research the products formed by the interaction between plants and insects. The optimum selection of stable RGs under temperature stress is necessary for follow up qRT-PCR experiments. In this study, we screened the expression of 11 candidate RGs: (Ubiquitin (UBQ), Tubulin1 (TUB1), TUB2, Tap interacting protein of 41 kDa (TIP41), proteins phosphatase 2A subnit A1 (PP2A1), PP2A2, glyceraldehyde-3-phosphate dehydrogenase 1 (GAPDH1), GAPDH2, actin (ACT), act-in1 (ACT1), and E3 ubiquitin-protein ligase UPL7 (UPL7) [11,12] for stability under temperature stress.

Materials
RCM saplings as the materials are from the Gallnut Breeding Base in Cili County in China. Thirty saplings were transplanted from the field to pots and cultured at 24°C and 60% humidity for 30 days. After growing more than 5 true leaves, 9 saplings with the same growth trend were selected and subjected to stress treatment using an artificial climate box. The plants were treated with the following stresses: (a) cold treatment, placed in a light incubator at 4°C long-days (16 h light, 8 h darkness); (b) heat stress, placed in a light incubator at 34°C for long-days (16 h of light, 8 h of darkness); and (c) the control group was grown for long-days (16 h light, 8 h darkness) in a culture chamber at 24°C. Leaves were harvested from control plants and plants treated under each stress condition for 3, 6 or 24 h. The leaves of three saplings from each treatment were blended as biological replicates, and three samples were taken altogether (Fig. S3). The samples were quickly frozen in liquid nitrogen and stored at À80°C for RNA extraction and subsequent experiments.

Total RNA extraction and cDNA synthesis
Total RNA was extracted from all samples using the Tiangen RNA prep Pure Plant Plus Kit plant (Tiangen, China) according to the manufacturer's instructions, The RNA concentration was measured using a NanoDrop 2000, and the RNA integrity was assayed by agarose gel electrophoresis (Fig. S4). Only RNA samples with an OD260/280 ratio between 1.8 and 2.2 and an OD260/230 ratio greater than 2.0 that showed three discrete bands of 28S, 18S, and 5S were used for cDNA synthesis. Synthesis of cDNA was performed using the NovoScript Plus All-in-one First-Strand cDNA Synthesis SuperMix (Novoprotein, China). According to the manufacturer's instructions, genomic DNA was directly removed by adding gDNA Purge. The reverse transcription system of each sample was as follows: RNA template (1 lg), gDNA Purge (1 lL), Supermix (10 lL) and RNase Free Water (up to 20 lL) at 50°C for 30 min and 75°C for 5 min.

Primer design and qPCR
The primer sequences for ACT1 were the same as those used previously [8]; the primers for all other reference genes (RGs) were designed according to transcriptome sequencing data from our group (unpublished data). First, based on the transcriptome data, common RGs such as GAPDH were searched, and then, the genes with the highest read count and FPKM (the highest expression) and lowest absolute value of the log2FC (the highest stability) were compared with Arabidopsis thaliana to determine the homology and identified as RCM reference bases. The sequences of these RGs are shown in Table S1. The annotations for these genes were compared with those in the TAIR Arabidopsis thaliana and NCBI databases. Specific primers were designed by SnapGene, and the specificity was analysed by NCBI primer blast. The primers were synthesized by Qingke (Beijing, China). The details of the genes and primers are shown in Table 1. The specificity of each primer pair was evaluated by amplification and dissolution curve analysis. The correlation coefficient (R 2 ) and amplification efficiency (E) of the primer pairs were evaluated using standard curves for mixed cDNAs that were diluted fivefold (1/5, 1/25, 1/125, 1/625 and 1/3,125).
The qPCR instrument (Roche LightCycler 480) used 384 well blocks. Each 20 lL reaction included 10.0 lL 29NovoStart SYBR qPCR SuperMix Plus (Novoprotein, China), 0.6 lL of each primer (10 lM), 2.0 lL diluted (1 : 10) cDNA template, and 6.8 lL RNase free water. The programme for product amplification was 94°C for 3 min followed by 45 cycles of 94°C for 15 s and 60°C for 35 s. The melting curve was generated after 45 cycles. Three technical replicates were performed for each sample. The steps of the melting curve were set at 95°C for 5 s (ramp rate: 4.80) and 65°C for 1 min (ramp rate: 2.50).

Data analysis
The amplification efficiency was calculated by E = 10 (À1/k) -1, where K is the slope of the standard curve [13]. The 2 ÀΔΔCt method was applied to calculate the relative expression level of the target gene, and significant difference analysis was conducted with the SPSS statistical software v23.0.

Algorithm to perform
geNorm was used to screen RGs, and eventually, 2 or more combinations of RGs were selected from the results The data input required by NormFinder is the same as that required by geNorm [15], which was processed to obtain the 2 ÀΔCp value for the RGs. The stability value of candidate genes was obtained, and it was determined that a candidate gene with a smaller stability value was more suitable to be used as an RG.
In BestKeeper, the Cp value is directly input into the algorithm, and the correlation coefficient (R), standard deviation (SD), and coefficient of variation (CV) are obtained after the algorithm is used [16].
The stability of gene expression was evaluated with geNorm, NormFinder, and BestKeeper.

Specificity evaluation and PCR efficiency analysis of primer amplification
In this study, 11 candidate RGs were identified based on RCM transcriptome data. The candidates were screened by analysing the genome data and gene annotation of the related species Camelina sativa, and alignment analysis of the Arabidopsis thaliana genes from TAIR (http://www.arabidopsis.org). Table 1 lists information for the genes, including the correlation description, primer sequence, gene number of the Arabidopsis homologue, amplicon length, PCR amplification efficiency, and correlation coefficient (R 2 ). The amplicons generated by qRT-PCR ranged from 92-200 bp. The specificity of the primers was confirmed by the presence of a single peak in the dissolution curve. The qRT-PCR dissolution curve peaks for the 11 RGs are shown in Fig. S1. Fig. 1 shows the crossing point (Cp) values for all samples (n = 33). The candidate gene with the highest abundance among all RGs was GAPDH2, which had the lowest mean Cp value of 13.34 AE 1.55 (mean AE SD) and thus the highest expression. In contrast, UPL7 (20.7 AE 0.971) had the highest mean Cp value (lowest expression). In addition, PP2A2 showed low variability and the lowest SD value (SD = 0.567) and might therefore be the most stably expressed RG in all samples, but further algorithms are required to support this hypothesis. Although the expression level of GAPDH2 (SD = 1.55) showed the most variability, and the SD was greater than 1, its unstable expression as an RG could reduce the accuracy of the results. Fig. S2 provides the separate Cp results of these genes for the control, heat, and cold treatments with respect to the time points (3, 6 and 24 h). To further evaluate the stability of the expression of the 11 genes in different samples exposed to different stresses, different algorithms were used to analyse the data.

Expression stability analysis of candidate reference genes
To further evaluate the stability of the expression of different RGs in response to different temperature stresses using other tools, the qRT-PCR data Cp values were corrected and preprocessed into appropriate

geNorm analysis
The whole dataset was divided into the following four groups: all samples group, cold stress sample group, heat stress sample group, and untreated control sample group. The results showed that the stability of the 11 RGs differed among these groups (Fig. 2). For example, GAPDH2 and UPL7 showed poor stability under various treatments and controls. ACT, ACT1, and PP2A2 were more stably expressed than other RGs in the control group (Fig. 2B) and in the heat stressed group (Fig. 2D). In the cold treated group, TUB2, GAPDH1, and PP2A1 showed highly stable expression (Fig. 2C). According to a Vn/Vn + 1 cutoff value of 0.15 for pairwise variations, the first number of n RGs is sufficient for accurate normalization [16]. In this article, the V2/3 values of the samples under cold and heat stress and the control group were less than 0.15 (Fig. 3), which indicates that the two RGs are sufficient for accurate normalization of all the samples. The Vn/n + 1 values of all samples were less than 0.15, except for V2/3 of the all-sample group and V10/11 of the heat stress group (Table S3), indicating that accurate normalization and reliable results can be obtained by using the n value corresponding to the number of RGs.

NormFinder analysis
Based on the results of the NormFinder analysis, Table 2 shows the expression stability and ranking of  (Table 3). In the heat-treated group, UBQ (1.1 AE 0.15), PIP2 (1.41 AE 0.28), and TIP41 (1.57 AE 0.31) were the top RGs. In addition, GAPDH1 and TIP41 were highly stable in the cold and heat-treated groups. In contrast, the SD value of GAPDH2 was greater than 1 in the heat-treated group, and the stability ranking in the other groups was also low, indicating that its stability was clearly the worst. The BestKeeper analysis suggested that the unstable expression of an RG is indicated by SD greater than 1 [16]. Although only two RGs, TUB1 under cold stress and GAPDH2 under heat stress, showed SD values greater than 1 in our analysis, half of the 11 candidate genes had CV values lower than 4.

Comprehensive stability ranking
The geometric means of different RG rankings are used to determine the most stably expressed RGs   Table 4. GAPDH1 (3.00), PP2A2 (3.33), and UBQ (3.33) ranked the top three in the heat stress group, while GAPDH1(2.33) and ACT (3.00) ranked among the top two in the cold stress group. The results from all the treatment groups indicate that GAPDH1 (3.00) and ACT (3.00) were the most reliable RGs among the 11 candidate genes. In the control (no-treatment) group, TIP41 (2.67) and GAPDH1 (2.33) were the most reliable RGs among the 11 candidate genes. GAPDH2 showed poor stability in all the treatment groups and was therefore unsuitable for use as an RG. It has good stability of GAPDH1 that the ranking in each result was ≤ 3.0. So GAPDH1 + PP2A2/UBQ were recommended for use in the heat stress group, and GAPDH1 + ACT were recommended for use in the cold stress group.

Validation of the selected reference genes
To validate the stability of the potential RGs selected using different algorithms and evaluation methods, the relative expression levels of OST1-1 and OST1-2 under temperature stress were verified, these are temperaturerelated genes of interest in the RCM transcriptome [19], and changes in the stability of candidate RGs under 3, 6, and 24 h of cold and heat stress were analysed based on their relative expression levels (Fig. 4). Under cold stress, ACT and GAPDH1 ranked as the most stable in the comprehensive evaluation and were used to verify OST after normalization by the two RGs. In Fig. 4A, the relative expression of OST1-1 was upregulated at 6 h compared with 3 h and was significantly upregulated at 24 h by approximately 11-fold. The expression of OST1-2 was upregulated at 6 h compared with 3 h and was significantly upregulated at 24 h by approximately 13-fold (Fig. 4B). Under heat stress, GAPDH1 and UBQ were selected as RGs with stable expression for verification. Although the average comprehensive ranking of UBQ and PP2A2 was the same under heat stress, the highest ranking of UBQ in the three algorithms was higher than that of PP2A2, so two pairs of RGs, UBQ and GAPDH1, were selected for verification. The expression levels of OST1-1 were upregulated by approximately 40 times and 80 times at 6 h and 24 h (Fig. 4C), respectively. The expression levels of OST1-2 were upregulated approximately 4-and 6-fold at 6 h and 24 h (Fig. 4D), respectively. When GAPDH2 with the lowest stability ranking was used as an RG to verify the relative expression level of the OST, the results obtained were completely opposite to the relative expression level of other candidate genes as an RG at 6 h, which was precisely caused by the instability of GAPDH2 itself. This result confirms the importance of the stability of RGs for experimental expression measurements under different conditions.

Discussion
Although the RCM genome has not been sequenced, RCM plants possess important medicinal and industrial value, which is mainly due to their high tannin contents [7,20,21]. Environmental temperature is a critical factor for RCM growth and gallnut formation.
Gallnuts are usually collected from RCMs planted on the shady side of mountains. The main goal of identifying stably expressed RGs in response to temperature in RCM is to reveal the basis of gallnut formation in  The stability of internal RGs according to geNorm analysis is indicated by the size of the M value; however, this value does not give an accurate threshold above which the internal RGs are unstable. The most stable RGs can only be selected as those with the lowest M value. Here, the M value for 11 candidate RGs was less than 0.8 under cold stress and less than 1.1 under heat stress, but the difference between the most stable and most unstable M value was 0.80 for heat stress treatments and only 0.39 for cold stress treatments. This suggests that the M value of the most unstable RG under cold stress was twice as high as that of the most stable RG, whereas that under heat stress was approximately 4-fold higher. A similar result was obtained from analysis using NormFinder, which also did not provide an accurate threshold with which to estimate the cut-off between candidate gene stability and instability.
In our study, among the 11 candidate RGs, 4 pairs of RGs (ACT1/ACT2, GAPDH1/GAPDH2, PP2A1/ PP2A2, and TUB1/TUB2) belonged to the same gene families but showed great differences in stability. As an example, the two candidate genes with the same name as GAPDH showed relatively large stability differences, which was also found in other studies [26]. Whether these differences are caused by the genes themselves or the differences caused by the specific amplified fragments that affect the amplification experiment is worth further study.

Conclusions
In this study, the recommended RGs under heat and cold stresses were similar to those recommended for 19 plant species. It is generally recommended to use multiple RGs to normalize expression data for most plant materials to improve accuracy. Because it is difficult to select appropriate RGs for each individual study, the statistics for different RGs here provide a good resource for gene expression analysis. Because differences in the expression of key genes are subtle, it is extremely important that RG expression stably responds to treatment. Therefore, this research provides information that is useful for the verification of gene expression in RCM under temperature stress and is useful for future research on gallnut formation. The screening and evaluation of the stability of the expression of 11 candidate RGs in RCM under temperature stress identified GAPDH1+ PP2A2/UBQ as the most suitable RG under heat stress, while GAPDH1+ACT are the most suitable under cold stress. In the untreated group, GAPDH1+TIP41, which are RGs, are recommended in general. In all treatment groups, GAPDH1+ACT were the recommended RGs. The results provide key technical guidance for the further study of RCMs under different temperature stresses and may also suggest candidate RGs for studies with other plant species.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article.    Data are mean AE SD and analysed by one-way ANOVA. Different letters above the bars represent significant differences (P < 0.05). Table S1. Sequence of RGs.