Identification and evaluation of reference genes for quantitative real-time PCR analysis in Polygonum cuspidatum based on transcriptome data

Background Polygonum cuspidatum of the Polygonaceae family is a traditional medicinal plant with many bioactive compounds that play important roles in human health and stress responses. Research has attempted to identify biosynthesis genes and metabolic pathways in this species, and quantitative real-time PCR (RT-qPCR) has commonly been used to detect gene expression because of its speed, sensitivity, and specificity. However, no P. cuspidatum reference genes have been identified, which hinders gene expression studies. Here, we aimed to identify suitable reference genes for accurate and reliable normalization of P. cuspidatum RT-qPCR data. Results Twelve candidate reference genes, including nine common (ACT, TUA, TUB, GAPDH, EF-1γ, UBQ, UBC, 60SrRNA, and eIF6A) and three novel (SKD1, YLS8, and NDUFA13), were analyzed in different tissues (root, stem, and leaf) without treatment and in leaves under abiotic stresses (salt, ultraviolet [UV], cold, heat, and drought) and hormone stimuli (abscisic acid [ABA], ethylene [ETH], gibberellin [GA3], methyl jasmonate [MeJA], and salicylic acid [SA]). Expression stability in 65 samples was calculated using the △CT method, geNorm, NormFinder, BestKeeper, and RefFinder. Two reference genes (NDUFA13 and EF-1γ) were sufficient to normalize gene expression across all sample sets. They were also the two most stable genes for abiotic stresses and different tissues, whereas NDUFA13 and SKD1 were the top two choices for hormone stimuli. Considering individual experimental sets, GAPDH was the top-ranked gene under ABA, ETH, and GA3 treatments, while 60SrRNA showed good stability under MeJA and cold treatments. ACT, UBC, and TUB were suitable genes for drought, UV, and ABA treatments, respectively. TUA was not suitable because of its considerable variation in expression under different conditions. The expression patterns of PcPAL, PcSTS, and PcMYB4 under UV and SA treatments and in different tissues normalized by stable and unstable reference genes demonstrated the suitability of the optimal reference genes. Conclusions We propose NDUFA13 and EF-1γ as reference genes to normalize P. cuspidatum expression data. To our knowledge, this is the first systematic study of reference genes in P. cuspidatum which could help advance molecular biology research in P. cuspidatum and allied species.


Background
Polygonum cuspidatum is a perennial herb of the Polygonaceae family that is a well-known traditional Chinese medicine. Its first medicinal application records date back over 1800 years in Mingyi Bielu [1]. It has long been used in Chinese folk medicine for the treatment of abdominal masses, postpartum blood stasis, urethritis, suppuration, ulcers, sore throats, toothache, chronic bronchitis, hemorrhoids, and other ailments [1,2]. More recently, oral intake or external applications of its processed products were shown to be effective in hepatitis, hypertension, hyperlipidemia, diabetes, jaundice, arthritis, skin burns, and scalds [3].
Numerous active compounds have been isolated and identified from P. cuspidatum to date, such as anthraquinones, flavonoids, stilbenes, organic acids, coumarins, catechins, and lignans [3]. The content of resveratrol, a stilbene, is higher in P. cuspidatum than in other plants [4], and several physiological and pharmacological studies have reported anti-oxidative, anti-cancer, antiinflammatory, anti-tumor, anti-depressant, and anti-viral roles for resveratrol, as well as neuroprotective and metabolic regulation through interactions with multiple targets [5]. Other components of P. cuspidatum also show significant health-promoting effects in cells or animal models [6,7]. Interestingly, the vast majority of active compounds are secondary metabolites synthesized during normal plant growth or in response to environmental stresses [8,9]. However, the mechanism of P. cuspidatum responses to abiotic stresses and hormone stimuli as well as in various tissues needs to be thoroughly explored to better understand the production of these active ingredients.
Stilbenes and flavonoids, derivatives of the phenylpropanoid pathway, actively participate in the regulation of resistance to stresses, including pathogens, ultraviolet (UV) radiation, low/high temperature, drought, heavy metals, methyl jasmonate (MeJA), and ethylene (ETH) [8,9]. Phenylalanine ammonia lyase (PAL) catalyzes the conversion of L-phenylalanine into cinnamic acid, which is the first committed step of the phenylpropanoid pathway. The same substrates 4-coumaroyl-CoA and malonyl-CoA are used by stilbene synthase (STS) and chalcone synthase to produce resveratrol and tetrahydroxychalcone, respectively, which are branch sites in the stilbene and flavonoid pathway [10]. Various transcription factors are involved in regulation of the phenylpropanoid pathway [11], and MYB4 in Arabidopsis responds to UV and SA stresses to suppress phenylpropanoid pathway gene expression [12,13].
The analysis of gene expression patterns can provide insights into complex metabolic processes. Quantitative real-time PCR (RT-qPCR) is commonly used to detect gene expression in different species because of its simplicity, speed, sensitivity, specificity, and high throughput [14]. However, despite these advantages, variability in initial materials, RNA integrity, RT-PCR efficiency, qPCR efficiency, and inherent technical variations will inevitably lead to errors, so it is necessary to use reference genes as internal controls [15]. An appropriate reference gene should be constantly expressed across the samples being investigated [15]. Classic reference genes for RT-qPCR data normalization, such as actin, tubulin, ubiquitin, elongation factor, translation initiation factor, and ribosomal RNA, are usually required for basic and essential processes in the cell. However, increasing experimental evidence shows that some of these genes are not as stable as previously thought [16]. Therefore, a host of highly stable novel reference genes have been identified using screening from genomes and transcriptome datasets [17,18]. Previously, reference genes have not always been stable in different species and tissues, or at different developmental stages and experimental conditions within a single species [19]. Therefore, it is crucial to assess expression stability probabilities under specific conditions prior to use. Furthermore, two or more reference genes are desirable to avoid a biased or mistaken interpretation of the changes of target gene expression [20].
In this context, the △CT method [21], geNorm [22], NormFinder [23], BestKeeper [24], and RefFinder [25] have been developed to evaluate the stability of genes in biological samples. They have been successfully employed to validate reference genes in various medicinal plants, such as Gentiana macrophylla [26], Achyranthes bidentata [27], and Euscaphis konishii [18]. In P. cuspidatum, PcPKS1 [28], PcPKS2 [29], PcCHS1 [30], PcSTS [31], PcMYB1 [32], and PcWRKY33 [33] have been cloned and studied. However, the lack of P. cuspidatum reference genes has become a major hurdle for gene expression studies in this species. Because of a lack of sequence and expression information, we performed transcriptome sequencing of P. cuspidatum vegetative tissues (root, stem, and leaf) and leaves under different treatments (UV and MeJA) in our laboratory (unpublished data), providing a wealth of resources for our current selection of suitable reference genes.
The present study was undertaken to characterize candidate reference genes for RT-qPCR data normalization in P. cuspidatum. ). Expression stability was calculated using △CT, geNorm, Norm-Finder, BestKeeper, and RefFinder. Additionally, the expression of three target genes (PcPAL, PcSTS, and PcMYB4) under UV and SA treatments and in different tissues was normalized separately by the most and least stable genes to demonstrate the suitability of the recommended reference genes. Our results provide some useful resources for the future quantification of gene expression in P. cuspidatum and allied species.

Identification of candidate reference genes
In this work, nine classic reference genes (ACT, TUA, TUB, GAPDH, EF-1γ, UBQ, UBC, 60SrRNA, and eIF6A) and three novel reference genes (SKD1, YLS8, and NDUFA13) were identified from the P. cuspidatum transcriptome as candidate reference genes. The full-length cDNAs of these 12 genes, which were used to design specific qPCR primers, were submitted to NCBI Gen-Bank ( Table 1).

Evaluation of primer specificity and efficiency
PCR products amplifying P. cuspidatum leaf cDNA were checked by 1% (m/v) agarose gel electrophoresis, which revealed a single band of expected size for each primer pair (Additional file 1: Figure S1). Additionally, a single peak in the melting curve for each gene confirmed primer specificity (Additional file 1: Figure S2). The qPCR efficiency (E) varied from 86% (TUB) to 103% (SKD1) with correlation coefficients (R 2 ) ranging from 0.994 to 0.999 (Additional file 1: Figure S3). Hence each primer pair was highly efficient and specific to the targeted region.

Expression profiling of candidate reference genes
The expression levels of the 12 candidate reference genes (ACT, TUA, TUB, GAPDH, EF-1γ, UBQ, UBC, 60SrRNA, eIF6A, SKD1, YLS8, and NDUFA13) were evaluated in 65 samples collected from different tissues without treatments and leaves under abiotic and hormone stimuli using CT values. Raw CT values of the 12 genes across all samples were shown in Additional file 2, and were plotted directly by Boxplot (Fig. 1). The CT values exhibited a relatively wide range from 19.84 (GAPDH) to 33.46 (SKD1) across all samples. Because the gene expression level is negatively correlated with the CT value, ACT was the highest expressed gene with the lowest mean CT value (21.60), whereas SKD1 was the least abundant gene with the highest mean CT value (27.42) among the 12 genes. The expression variation among the 65 samples for each gene ranged from 4.12 (eIF6A) to 8.779 (SKD1). No gene showed an unchanged expression level under all conditions, so it was necessary to identify reference genes under specific experimental conditions in P. cuspidatum.

Expression stability of candidate reference genes
In our study, each reference gene was evaluated in 11 experimental sets which were analyzed individually at first. To obtain a more comprehensive analysis, the sets were then divided into four different groups: (1) "Abiotic stress" (salt, UV, cold, heat, and drought) (2) "Hormone stimuli" (ABA, ETH, GA 3 , MeJA, and SA) (3) "Different tissues" (root, stem, and leaf), and (4) "All" (all experimental sets). More specifically, the stability of the 12 candidate reference genes was analyzed by the △CT method, geNorm, NormFinder, BestKeeper, and RefFinder.
geNorm estimates the optimal number of reference genes required for accurate and reliable RT-qPCR normalization. As shown in Fig. 2, V 2/3 values were much lower than the cut-off value of 0.15 for all different groups, indicating that two reference genes would be sufficient for accurate and reliable normalization of the gene expression data.
For "Abiotic stress", 60SrRNA and EF-1γ had the lowest M value (0.25) followed by NDUFA13 (0.35), indicating that these were the three most stable genes in geNorm analysis. NDUFA13 and EF-1γ were the two top-ranking genes in NormFinder and △CT analysis, while eIF6A and 60SrRNA were top in BestKeeper analysis. According to the calculations performed by RefFinder, EF-1γ and NDUFA13 were the most appropriate reference genes under abiotic stress conditions. The comprehensive ranking order for every specific abiotic stress was not entirely consistent. EF-1γ, NDUFA13, and 60SrRNA were also identified as the three most stable reference genes under heat treatment, but EF-1γ and NDUFA13 fell outside the list of the four most stable genes, though 60SrRNA ranked top under cold conditions. For salt treatment, YLS8 was the top-ranked gene, followed by EF-1γ, TUB, and NDUFA13. For UV treatment, UBC had the highest stability, while EF-1γ, NDUFA13, and 60SrRNA also exhibited good stability. NDUFA13 was the optimal reference gene under drought treatment, followed by EF-1γ and 60SrRNA with low CVs (Table 2, Additional file 3: Table S1).
For "Hormone stimuli", the best reference genes were 60SrRNA and EF-1γ in geNorm analysis. NDUFA13 and SKD1 were the two most stable reference genes by both △CT and NormFinder, but ACT and SKD1 were identified by BestKeeper. In a comprehensive analysis, NDUFA13 and SKD1 were the two optimal reference genes under hormone stimuli conditions. When considering every hormone treatment, NDUFA13 was one of the three most stable reference genes for all hormone treatments except ABA; GAPDH and SKD1 were among the four most stable genes under ABA, ETH, and GA conditions ( Table 2, Additional file 3: Table S1).
For "Different tissues", NDUFA13 and EF-1γ were the most stable combination in geNorm analysis. Similar results were seen by △CT. NDUFA13 and EF-1γ next to GAPDH were the top-ranking genes in NormFinder analysis. In BestKeeper analysis, ACT ranked first with the lowest stability value (0.06), while NDUFA13 and EF-1γ ranked second and fourth, respectively. NDUFA13 and EF-1γ were recommended in different tissues of P. cuspidatum by RefFinder (Table 2, Additional file 3: Table  S1).
When all samples were taken into account, 60SrRNA, EF-1γ, and NDUFA13 ranked most highly in geNorm, △CT, and NormFinder analysis. However, they were low in the ranking in BestKeeper analysis. RefFinder ranked the candidate reference genes from the highest to the lowest stability as follows: Table 2). Taken together,  Table S1).

Validation of candidate reference genes
To ensure the accuracy and reliability of our results, the relative expression patterns of PcPAL, PcSTS, and PcMYB4 were analyzed under Abiotic stress (UV), Hormone stimuli (SA), and in different tissues (leaf, stem, and root). The two most stable reference genes (NDUFA13 and UBQ for SA, UBC and EF-1γ for UV, NDUFA13 and EF-1γ for different tissues) and one unstable gene (TUA) were selected for normalizing qPCR data. As shown in Fig. 3a, under UV treatment, PcPAL was significantly induced at all analyzed stress times. However, the expression levels at 16 h and 32 h with TUA normalization were much higher than those with either UBC, EF-1γ, or their geometric mean. The expression of PcMYB4 was inhibited by UV treatment, initially reduced after reaching the lowest level at 8 h, then increased afterwards. However, the expression levels of PcMYB4 at 16 h and 32 h were overestimated with TUA normalization. Similar misjudgments occurred in the analysis of PcSTS expression data. Some divergences in the results were also   observed under SA treatment (Fig. 3b). The expression level of PcPAL with NDUFA13, UBQ, or their geometric mean achieved the highest after adding SA for 2 h and was higher at 4 h than the control, whereas the opposite results were displayed with TUA normalization. PcSTS and PcMYB4 responded quickly to SA, and a low expression level was maintained for a long time (2-12 h). TUA normalization similarly led to erroneous interpretation of the relative expression patterns of PcMYB4. In different tissues (Fig. 3c), the expression patterns of PcPAL, PcSTS, and PcMYB4 were similar when NDUFA13,  The most stable reference genes (NDUFA13 and UBQ for SA, UBC and EF-1γ for UV, NDUFA13 and EF-1γ for tissues) and unstable gene (TUA) were selected for normalizing qPCR data. Data are represented as mean ± SD, different letters on the vertical bars indicate significant difference at 0.05 levels EF-1γ, and TUA were used for normalization, but the fold-changes in root and stem were underestimated with normalization by TUA.
Overall, the expression patterns of PcPAL, PcSTS, or PcMYB4 were nearly the same when using stable reference genes for normalization, whereas the expression levels had large variations when TUA was used. Furthermore, the RT-qPCR results normalized by the stable reference genes under UV treatment and in different tissues were more consistent with the target gene expression profiling derived from P. cuspidatum transcriptome data (Additional file 3: Table S2). Thus, the most suitable reference genes calculated by the above mentioned software were applicable under some specific conditions.

Discussion
RT-qPCR is one of the most commonly used techniques to obtain gene expression profiles in molecular biology. A prerequisite of this is the selection of appropriate reference genes for data normalization to ensure the accuracy and reliability of the results [15]. Large-scale gene segments and gene expression data generated by sequencing provide abundant resources for the identification and evaluation of reference genes, especially in nonmodel species [18,26].
In this work, we made full use of transcriptome sequencing data available in our laboratory to identify candidate reference genes for P. cuspidatum. The expression stability was evaluated by △CT method, geNorm, NormFinder, and BestKeeper, although the results obtained were not completely consistent. Such discrepancies in stability ranking were also reported in some previous studies [34]. Interestingly, we found that the rankings by △CT, geNorm, and NormFinder, were similar, especially for individual sets, but different from those by BestKeeper; similar findings were reported in other studies [19,34]. For example, in our work, under ETH, GA 3 , and SA treatment, ACT was identified as the most stable gene by BestKeeper analysis, but performed unsatisfactorily in △CT, NormFinder, and geNorm analysis. Therefore, comprehensive analysis with RefFinder of multiple results from different software could help select a more appropriate reference gene.
ACT, TUA, and TUB, which encode cytoskeletal proteins, are extensively used as reference genes, and their high stability is reported in many plants [35,36]. However, we showed that TUA performed poorly across all sample sets, ACT performed less well in most cases, and TUB was more affected by environmental factors such as SA and cold; this was similar to previous results obtained in Eucommia ulmoides [37]. GAPDH encodes an abundant glycolytic enzyme present in most cell types, and has universally been used as a reference gene in RT-qPCR. GAPDH showed a good performance in response to "Hormone stimuli" in parsley and carrot leaves, but was less stable for most abiotic stresses [36,38]. The performance of GAPDH in our study was relatively consistent with those of previous studies, except for responses to MeJA and drought. Ubiquitin is a small regulatory protein found in most tissues of eukaryotic organisms. We showed that UBC was the most stable gene under UV treatment, in agreement with the findings of Borges [39]. Moreover, UBQ was one of the two most stable genes under SA treatment, which differed from the results obtained in Achyranthes bidentata [27] and carrot [36]. However, UBC and UBQ did not perform well across all other sets, which is consistent with the results seen in parsley [38]. Thus, UBC and UBQ were not the preferred gene choices in our study.
Eukaryotic translation initiation factor and ribosomal RNA have been reported as reference genes in many studies, such as eIF4a in Eleusine indica [40], and eIF5A and 60S ribosomal RNA in Panax ginseng [41]. In our study, the performances of eIF6A and 60SrRNA were not always the best in "Abiotic stress", "Hormone stimuli", and "all" groups, but were within an acceptable range. Elongation factor 1 (EF-1), composed of the four subunits EF-1α, EF-1β, EF-1δ, and EF-1γ, plays a central role in protein biosynthesis [42]. EF-1α, encoding a Gprotein, was used for the normalization of qPCR data in some medicinal plants such as G. macrophylla [26] and A. bidentata [27]. The main function of EF-1γ is to ensure the correct scaffolding of different subunits in the EF-1 complex as well as to direct its intracellular localization [43]. The expression stability of EF-1γ had previously been evaluated in P. ginseng [34] and Nilaparvata lugens [44]. In our study, EF-1γ was one of the two most suitable reference genes with stable expression levels under various experimental conditions.
Novel reference genes such as SKD1, YLS8, and NDUFA13 have also been selected for gene normalization by the application of large amounts of omics data. SKD1, encoding a protein that contributes to vacuolar trafficking and maintenance of the large central vacuole of plant cells [45], was first selected as a reference gene in pear [46]. In this study, SKD1 was preferable under "Hormone stimuli" conditions, but had to be discarded under "Abiotic stress" and "Different tissues" conditions because of its low expression and high variation. YLS8 encodes a protein involved in mitosis, and performed well under ABA, MeJA, and salt conditions, but poorly in other individual sets and the three combination groups. These results were not identical to those obtained previously [46,47]. Notably, NDUFA13 was remarkably stable in the individual sets (except in the ABA group) as well as in the combination groups, suggesting it is an almost ideal reference gene. NDUFA13, encoded by the nuclear genome, is an accessory subunit of the mitochondria respiratory chain complex I, which transfers electrons from nicotinamide adenine dinucleotide to ubiquinone [48]. The fact that mitochondria are known as eukaryotic cell powerhouses may explain why NDUFA13 showed stable expression under various conditions in this study. NDUFA13 was previously reported as a reference gene in Apostichopus japonicus [49], but no studies have investigated its potential in plants. In this regard, the present study is the first known report of NDUFA13 as a reference gene in plants.
PAL, STS, and MYB4 are three important regulating genes in the phenylpropanoid pathway. When using the stable reference genes, the relative expression patterns of PcPAL and PcMYB4 were in agreement with previous reports; for example, MYB4 was down-regulated by exposure to UV-B light in Arabidopsis thaliana [12] and Brassica rapa [50], and PAL was induced by exposure to UV and SA in A. thaliana [51] and Juglans regia [52]. To our knowledge, STS is found in only a few plants, and not in Arabidopsis or tobacco. Research into STS expression is currently concentrated in grapes, where it was shown to be strongly induced by UV-C [53] and SA [54]. We observed the opposite expression patterns of STS in response to SA in P. cuspidatum, with rapidly decreasing expression lasting for 12 h. These differences may reflect variations in UV wavelengths or species. Moreover, our results are consistent with P. cuspidatum transcriptome data. Evidently, there were obvious underestimates or overestimates in our results following normalization by the least stable gene, TUA. Our results demonstrate the importance of using a stable reference gene for normalization to obtain accurate results.
Based on the above analysis, we suggest that no single gene should be used for normalization in all species, tissues, or treatments. Therefore, suitable reference genes for given species and conditions should be explored.

Conclusions
We evaluated the expression stability of 12 candidate reference genes in different tissues of P. cuspidatum and under different treatment conditions. NDUFA13 and EF-1γ were identified as the best two reference genes for normalizing RT-qPCR gene expression data. Their reliability and effectiveness were verified by PcPAL, PcSTS, and PcMYB4. To our knowledge, the current work is the first systematic analysis of suitable reference genes that will facilitate further research into the molecular biology of P. cuspidatum and other closely related species.

Plant material
Seeds of P. cuspidatum were collected from the medicinal plant garden of the Institute of Botany, Chinese Academy of Sciences (Beijing, China). Seeds were surface-sterilized and sown on Murashige and Skoog (MS) agar medium in a growth chamber at 24°C with a 16 h/8 h light/dark cycle. After 1 month, similar seedlings were submitted to different treatments. For hormone treatment, the seedlings were uniformly sprayed with 0.5 mM MeJA, 1 mM SA, 0.1 mM ABA, 0.5 mM GA 3 or 0.5 g/L ethrel. For cold and heat treatments, the seedlings were respectively kept in 4°C and 42°C illumination boxes. For UV irradiation, the seedlings were placed under a UV-B transilluminator for 20 min. For salt treatment, the seedlings were transplanted to MS agar medium containing 100 mM NaCl. Leaf samples were collected at 0, 2, 4, 8, 16, and 32 h after UV and MeJA treatments, and 0, 2, 4, 8, 12, and 24 h for the other treatments. For drought treatment, plants grown for 2 months in plastic pots containing a soil/vermiculite mixture (1:1) under the same conditions were kept in dry soil for 1 week, then rehydrated and sampled daily (0, 1, 2, 3, 4, 5, 6, 7, and 8 d). Tissue-specific samples (root, stem, and leaf) were collected from plants grown for 2 months in the soil. All 65 samples, including three tissue-specific samples (root, stem, and leaf), 32 stresstreated samples (salt-, UV-, cold-, heat-, and droughttreated leaves), and 30 hormone-treated samples (ABA-, ETH-, GA 3 -, MeJA-, and SA-treated leaves), were separately collected in three biological repeats. All samples were immediately frozen in liquid nitrogen and stored at − 70°C before RNA extraction.

RNA isolation and cDNA synthesis
Total RNA was isolated from samples using the Plant Total RNA Purification Kit (GeneMark, Taiwan, China) following the manufacturer's instructions. The RNA integrity was checked on a 1% (m/v) agarose gel. The quantity and quality of the total RNA samples were assessed by recording absorbances at 260/280 nm and 260/230 nm with the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Only RNA samples 1.8 < OD 260/280 < 2.2 and OD 260/230 > 2.0 were used for subsequent cDNA synthesis. Total RNA (3 μg) was reverse-transcribed into cDNA using the Hifair™ II 1st Strand cDNA Synthesis SuperMix (Yeasen, Shanghai, China) with oligo (dT) primers according to the manufacturer's instructions. cDNA was diluted at 1:20 with the EASY dilution solution (Takara, Japan), then stored at − 20°C until required as template for qPCR.

Selection of candidate reference genes
The expression levels of unigenes are commonly estimated by fragments per kilobase of transcript per million mapped reads values [55]. A candidate reference gene should be moderately expressed with a small coefficient of variation (CV) [56]. In this work, the CV calculation formula was: CV = standard deviation of reads per kilobase of transcript per million mapped reads (RPKM)/average RPKM. We downloaded the sequences of nine classic reference genes (ACT, TUA, TUB, GAPDH, EF-1γ, UBQ, UBC, 60SrRNA, and eIF6A) identified in other species, and carried out BLASTn queries against the P. cuspidatum transcriptome. Then, the highest ortholog sequences were checked against the Arabidopsis genome database and recorded in Table 1. Additionally, three new candidate genes, SKD1, YLS8, and NDUFA13, were selected based on their low CV and appropriate RPKM values (Additional file 3: Table S3). The cDNA and genomic DNA sequences of these candidate reference genes were shown in Additional file 4 and Additional file 5, respectively.

Primer design
Specific primers were designed using Primer 5 software and synthesized by GENEWIZ Company (Tianjin, China). Gene characteristics and primer details are shown in Table 1. The specificity of each primer pair was assessed by amplification and melting curve analysis. The correlation coefficient (R 2 ) and amplification efficiency (E) for all primer pairs were evaluated by standard curves using fourfold dilutions of the pooled cDNA (1/4, 1/16, 1/64, 1/256, and 1/1024).

RT-qPCR
All qPCR was carried out in 96-well plates using Quant-Studio™ Real-Time PCR Software (Applied Biosystems, USA) with the Hieff™ qPCR SYBR® Green Master Mix (Yeasen, Shanghai, China). Each 10 μL reaction included: 5.0 μL of 2× Hieff™ qPCR SYBR® Green Master Mix, 0.2 μL of each primer (10 μM), 1.0 μL of diluted (1:20) cDNA template, and 3.6 μL of RNase-free water. The cycle program for product amplification was: 95°C for 5 min (hot-start activation) followed by 40 cycles of 95°C for 10 s (denaturation), 58°C for 20 s (annealing), and 72°C for 20 s (extension). The melting curve was generated after 40 cycles to test the specificity of each primer pair across the temperature range of 60-95°C at a heating rate of 0.05°C/s. Three technical replicates were used for each sample.
The △CT method was employed to rank the genes by calculating the average standard deviation (SD) based on the relative expression of all pairwise combinations of candidate reference genes. The gene with the lowest SD was identified as the most stable reference gene [21]. The geNorm calculates the expression stability value (M-value) for each gene. Genes with the lowest M values have the most stable expression. geNorm determines the pairwise variations (V) of one specific gene with all others. A cut-off value of V n/n + 1 < 0.15 means that further addition of reference genes no longer makes any significant contribution to the normalization. For instance, V 2/3 < 0.15 means that two reference genes are sufficient for data normalization [22]. NormFinder provides a stability value (SV) for each gene, which takes intergroup and intragroup relationships into consideration. A lower SV indicates a higher stability [23]. BestKeeper calculates three variables based on the CT values of all genes: the coefficient of correlation (r), standard deviation (SD), and coefficient of variance (CV). A more stable gene exhibits a lower SD ± CV value [24]. RefFinder gives a comprehensive ranking for each candidate reference gene based on the geometric mean of the weights of all genes calculated by the above four computational approaches. A lower geometric mean of the ranking values indicates a more stable expression [25].

Validation of selected candidate reference genes
To validate the selected reference genes, qPCR was performed to analyze target gene expression levels (PcMYB4, PcPAL, and PcSTS) under UV and SA treatments and in different tissues. The homologous P. cuspidatum genes AtMYB4 and AtPAL1 were cloned and named PcMYB4 and PcPAL, respectively (Table 1, Additional file 4, Additional file 5). PcSTS (EU647245.1) was isolated in previous studies from our laboratory [31]. Primer design and detection for these three genes were performed according to the aforementioned methods. Relative expression levels were calculated with the 2 −△△CT method.
Additional file 2: Raw CT values of the 12 candidate reference genes across all samples of Polygonum cuspidatum.