Identification and Validation of Reference Genes for RT-qPCR Analysis in Switchgrass under Heavy Metal Stresses

Switchgrass (Panicum Virgatum L.) has been recognized as the new energy plant, which makes it ideal for the development of phytoremediation on heavy metal contamination in soils with great potential. This study aimed to screen the best internal reference genes for the real-time quantitative PCR (RT-qPCR) in leaves and roots of switchgrass for investigating its response to various heavy metals, such as cadmium (Cd), lead (Pb), mercury (Hg), chromium (Cr), and arsenic (As). The stability of fourteen candidate reference genes was evaluated by BestKeeper, GeNorm, NormFinder, and RefFinder software. Our results identified U2AF as the best reference gene in Cd, Hg, Cr, and As treated leaves as well as in Hg, Pb, As, and Cr stressed root tissues. In Pb treated leaf tissues, 18S rRNA was demonstrated to be the best reference gene. CYP5 was determined to be the optimal reference gene in Cd treated root tissues. The least stable reference gene was identified to be CYP2 in all tested samples except for root tissues stressed by Pb. To further validate the initial screening results, we used the different sets of combinatory internal reference genes to analyze the expression of two metal transport associated genes (PvZIP4 and PvPDB8) in young leaves and roots of switchgrass. Our results demonstrated that the relative expression of the target genes consistently changed during the treatment when CYP5/UBQ1, U2AF/ACT12, eEF1a/U2AF, or 18S rRNA/ACT12 were combined as the internal reference genes. However, the time-dependent change pattern of the target genes was significantly altered when CYP2 was used as the internal reference gene. Therefore, the selection of the internal reference genes appropriate for specific experimental conditions is critical to ensure the accuracy and reliability of RT-qPCR. Our findings established a solid foundation to further study the gene regulatory network of switchgrass in response to heavy metal stress.


Introduction
Real-time quantitative PCR (RT-qPCR) has become the leading technique applied in gene expression analysis due to its advantageous characteristics, such as high-throughput, high sensitivity, and specificity, along with great repeatability. Reference genes with stable expression levels are used as the standard markers to calibrate and ensure the accuracy and validity of results from RT-qPCR [1,2]. Thus, a few conventional housekeeping genes, such as β-Actin and 18S rRNA, are commonly used as reference genes in RT-qPCR for plants [3]. However, based on up-to-date studies, the universal reference gene with stable expression profiles in different tissues and organs, developmental stages as well as experimental conditions has not been discovered. Therefore, it is of pivotal importance to identify the appropriate and stable reference genes associated with various situations in the analysis of gene expression profiles.
Heavy metal, as the main environmental pollutant, has raised a growing concern in ecological and global public health in recent years [4]. With the overdevelopment and excessive utilization of mineral resources and reserves, the disposal of electronic waste, and the extensive application of pesticides and fertilizers, the pollution in soil and water by heavy metals has become an increasingly serious problem [5]. Furthermore, the enrichment of heavy metals in the soil by agronomic crops endangers human health in a direct manner. Therefore, the remediation of heavy metal polluted soils, especially phytoremediation, has gained increasing attention from both academia and industries due to its lower cost and fewer side effects than conventional chemical and physical techniques [6]. Though the hyperaccumulator has become the hotspot for ecological restoration in recent years, utilization of energy plants for remedying heavy metal contaminated soil can achieve a win-win situation for the production of biomass raw materials as well as the management of the polluted soil [7].
Switchgrass (Panicum Virgatum L.), a C4 warm-season perennial grass species, is originated from North America and widely distributed in non-natural areas (above 55 • north latitude) [8]. The C4 grass has many merits, such as high efficiency in photosynthesis, high utilization of nitrogen, water, and nutrients, and it is also effective in water-soil conservation as well as the increase in organic matter in the soil [9]. Generally, the full growth season of switchgrass starts in the third year after plantation, and a single plantation lasts from 10 to 15 years. In addition, switchgrass is well tolerated and grows well on the land under a variety of stress conditions, such as drought, alkali-salt, flooding, and leanness [10]. The renewable energy that is produced from switchgrass has been reported to be about five fold of the energy that is consumed during the production process [11]. In addition, considering the environmental benefits in soil conservation and reducing greenhouse gas emissions, switchgrass has been recognized as the new energy plant with great potential, which makes it ideal for the development of phytoremediation on heavy metal contamination in soils [7,10].
In this study, we have selected fourteen common housekeeping genes as candidate reference genes from previous studies [10,12], aiming to identify appropriate reference genes with stable expression in various tissues of switchgrass and analyze its response to stress induced by different heavy metals. We further studied the differential gene expression of PvZIP4 and PvPDR8 from Zinc/Iron regulatory transporter family (ZIP) and ATP-binding cassette transporter family (ABC), respectively, pre and post the stress. Our results not only facilitate the understanding of the molecular mechanism of the heavy metal stress-induced response in switchgrass, but also establish the foundation for further studies on remediation of heavy metal contamination in soil.

Plant Materials and Treatment
The switchgrass cultivars "Alamo" seeds were sterilized in 10% NaClO for 30 min. After five times rinsing with deionized water, 100 seeds were sown and germinated in trays (30 cm length, 16 cm width, and 12 cm deep) with 1/2 strength Hoagland's solution. Plants were incubated in a growth room with the following environmental conditions: temperature 25 • C/20 • C (12 h day/12 h night), photosynthetically active radiation 400 µmol m −2 s −1 . Two months after plantation, plants were treated with 1 mM solution of cadmium (Cd), lead (Pb), mercury (Hg), chromium (Cr), and arsenic (As), respectively. Three replicates (three pots) were conducted for each heavy metal treatment group in a completely randomized design. Leaf and/or root samples were collected at the following time points 0, 1.5, 3, 6, 12, 24, and 48 h post-treatment with different heavy metals. Tissue samples were stored at −80 • C for further analysis.

Extraction of Total RNA and Reverse Transcription
Switchgrass tissues weighing 50-100 mg were pulverized using a Tissuelyzer (Qiagen, Germantown, MD, USA), and total RNA was extracted via Trizol reagent (Invitrogen, Carlsbad, CA, USA). The quality of RNA was validated by running the 1.2% agarose gel electrophoresis and visualized by an Analytikjena ScanDrop (Jena, Germany). Total RNA (1 µg) was extracted for each sample and reverse transcribed into the first strand of cDNA using the PrimeScript TM RT reagent kit (TaKaRa, Dalian, China) and stored at −80 • C for further analysis.

Real-Time PCR
Quantitative analysis via real-time PCR was conducted using a Sosofast Supermix reagent kit (Bio-Rad, Hercules, CA, USA) and the CFX96 real-time PCR system (Bio-Rad, Hercules, CA, USA). The experiment was performed in 20 µL volume of reaction in ice-bath. Reaction reagents were as follows: 1 µL of primer at a final concentration of 0.2 µmol·L −1 , 10 µL qPCR SYBR Green SuperMix, 2 µL cDNA, and ddH 2 O to bring the total volume to 20 µL. The sequential steps of real-time PCR include pre-denaturation at 94 • C for 10 s, denaturation at 94 • C for 10 s, annealing at 62 • C for 5 s, 40 cycles in total. Three technical replicates were performed in a sample mixture with each heavy metal stress at each time point. Then the two closest cycle threshold (Ct) values were used for RT-qPCR analysis.

Analysis of Stability
The cycle threshold (Ct) value for each reference gene was obtained from RT-qPCR and analyzed through GeNorm [13], NormFinder [14], BestKeeper [15], and RefFinder (http://www. leonxie.com/referencegene.php) software. When the data analysis was performed with GeNorm and NormFinder, the Ct value was first converted to relative quantitative Q value via the formula Q = 2 −∆Ct (∆Ct = Ct sample − Ct min ). Ct sample is the Ct value of the housekeeping gene in each of heavy metals treatment; Ct min indicates the lowest Ct value of this housekeeping gene among each of heavy metals treatment. Then the expression stability measurement (M) value was calculated by the GeNorm program for each candidate reference gene. BestKeeper directly utilized the Ct value for stability analysis without the additional converting step to measure the comparisons of the coefficient of variance (CV) and the standard deviation (SD). Finally, RefFinder integrates all three methods mentioned above to calculate the geometric mean for each reference gene and the comprehensive ranking index of stability. A lower index value indicates a higher stability of the reference gene. The optimal number of reference gene was determined by the paired coefficient of variation V n /V n+1 . It is generally considered that when the value of V n /V n+1 is less than 0.15, it is unnecessary to introduce a new reference gene. Otherwise, the (n+1) th reference gene is in need.

Validation of Reference Genes by Expression Analysis of Two Metal Transporters Genes
The two homologs of metal transporters PvZIP4 and PvPDR8 from switchgrass were obtained from the database (https://phytozome.jgi.doe.gov/pz/portal.html). For the validation of selected reference genes, the expression levels of two genes were analyzed using the most and least stable reference genes under heavy metal stress, calculated using the 2 −∆∆Ct method [16]. Three replicate samples were included for each treatment, and three technical replicates were conducted for each biological sample.

Specificity of Primers for Reference Genes
The RT-qPCR reactions were performed using the total RNA reverse transcription products from young leaves and young roots of switchgrass treated by different heavy metals treatment as the template. The results suggested that a distinct single peak was identified in the melting curves of all genes ( Figure S1). In addition, the PCR amplification curve of all samples showed great repeatability, indicating that the primers were able to amplify the desired products of each gene with high specificity and no primer dimer. Therefore, our results from RT-qPCR were confirmed to be valid and reliable.

Analysis of Reference Gene Expression
It has been reported that the Ct value of reference genes is inversely proportional to the expression level of that gene [17]. The greater the Ct value of the reference gene is, the lower the amount of the target gene being expressed in the sample and vice versa. Expression abundance of fourteen reference genes in all samples was analyzed via RT-qPCR ( Figure 1; Table S2). Our results demonstrated that the Ct value for each reference gene was in the range of 5-33. The lowest Ct value was found in the 18S rRNA gene ranging between 5 and 15, while its expression abundance was the highest among all the reference genes. The low values were in the case of UBQ2 and ACT12. ACT2 had high value as CYP2 and eEF4a, indicating the lowest expression abundance. The large distribution of Ct values suggested that the expression abundance differs among the reference genes. Median cycle threshold (Ct) values for fourteen candidate reference genes in switchgrass root and leaf samples under heavy metals stress conditions. The variation is displayed as medians values (lines across the box plot), 25th to 75th percentiles (boxes), and the maximum and minimum values (whiskers). Cadmium-treated leaves (CdL) and roots (CdR); lead-treated leaves (PbL) and roots (PbR); mercury-treated leaves (HgL) and roots (HgR); chromium-treated leaves (CrL) and roots (CrR); arsenic-treated leaves (AsL) and roots (AsR), the same below.

GeNorm Analysis
The expression stability of reference genes was analyzed via GeNorm V3.4 and represented by calculated M values. The lower the M value is, the higher stability the reference gene has and vice versa. Based on this principle, the M value was determined for each reference gene of all samples. Different combination of reference genes was shown to be the most stable ones in roots and leaves responding to each heavy metal stress: UBQ1/FTSH4 in Cd-treated roots (CdR) U2AF/ACT12 in Cd- Figure 1. Median cycle threshold (Ct) values for fourteen candidate reference genes in switchgrass root and leaf samples under heavy metals stress conditions. The variation is displayed as medians values (lines across the box plot), 25th to 75th percentiles (boxes), and the maximum and minimum values (whiskers). Cadmium-treated leaves (CdL) and roots (CdR); lead-treated leaves (PbL) and roots (PbR); mercury-treated leaves (HgL) and roots (HgR); chromium-treated leaves (CrL) and roots (CrR); arsenic-treated leaves (AsL) and roots (AsR), the same below.

GeNorm Analysis
The expression stability of reference genes was analyzed via GeNorm V3.4 and represented by calculated M values. The lower the M value is, the higher stability the reference gene has and vice versa. Based on this principle, the M value was determined for each reference gene of all samples. Different combination of reference genes was shown to be the most stable ones in roots and leaves responding to each heavy metal stress: UBQ1/FTSH4 in Cd-treated roots (CdR) U2AF/ACT12 in Cd-treated leaves (CdL), eEF1a/U2AF in Pb-treated roots (PbR), 18Sr/ACT12 in Pb-treated leaves (PbL), UBQ1/UCE in Hg-treated roots (HgR), UBQ6/CYP5 in Hg-treated leaves (HgL), UCE/UBC in Cr-treated roots (CrR), UCE/UBC in Cr-treated leaves (CrL), eEF1a/eEF4a in As-treated roots (AsR) and U2AF/ACT12 in As-treated leaves (AsL). However, the overall evaluation suggested that eEF4a and U2AF in both leaf and root tissues displayed the highest stability with the lowest M values under all stress conditions tested ( Figure 2).
Genes 2020, 11, x FOR PEER REVIEW 6 of 19 treated leaves (CdL), eEF1a/U2AF in Pb-treated roots (PbR), 18Sr/ACT12 in Pb-treated leaves (PbL), UBQ1/UCE in Hg-treated roots (HgR), UBQ6/CYP5 in Hg-treated leaves (HgL), UCE/UBC in Crtreated roots (CrR), UCE/UBC in Cr-treated leaves (CrL), eEF1a/eEF4a in As-treated roots (AsR) and U2AF/ACT12 in As-treated leaves (AsL). However, the overall evaluation suggested that eEF4a and U2AF in both leaf and root tissues displayed the highest stability with the lowest M values under all stress conditions tested ( Figure 2).  When the pairwise variation V n /V n+1 value is lower than the threshold of 0.15, the value (n) can be considered as the optimal number of reference genes for accurate normalization [13]. The V 2/3 value of reference gene in all samples under the stress of each heavy metals was shown to be smaller than the threshold value 0.15 ( Figure 3), indicating that the gene expression analysis needs two reference genes to achieve the best performance. However, the combined use of the four reference genes could be suitable for testing all the considered tissues and stress conditions. When the pairwise variation Vn/Vn+1 value is lower than the threshold of 0.15, the value (n) can be considered as the optimal number of reference genes for accurate normalization [13]. The V2/3 value of reference gene in all samples under the stress of each heavy metals was shown to be smaller than the threshold value 0.15 ( Figure 3), indicating that the gene expression analysis needs two reference genes to achieve the best performance. However, the combined use of the four reference genes could be suitable for testing all the considered tissues and stress conditions.

BestKeeper Analysis
The standard deviation (SD) of the Ct value of each housekeeping gene was calculated via BestKeeper V1. With the SD value less than 1, the housekeeping gene is considered as the stable one. Furthermore, the lower the SD value is, the higher stability that gene displays. Reversely, the gene was counted as being not stable if the SD value is higher than 1. Our results demonstrated that the U2AF gene exhibited the highest stability in total, particularly in the root sample treated with Hg and leaf samples treated with Pb, Hg, and As (Table 1). UCE was found to be the most stable reference gene in leaf samples treated with Cd, while eEF4a was identified to be the most stable reference gene for root samples treated with Cd and As. In the root and leaf samples treated with Pb and Cr, respectively, the UBC gene was shown to have the highest stability. When the root samples were treated with Cr, FTSH4 was found to be the most stable gene. Among all of the reference genes, the SD values of the CYP2 gene were demonstrated to be greater than 1 in all treatments except for PbR, indicating the instability of this gene (Table 1).

BestKeeper Analysis
The standard deviation (SD) of the Ct value of each housekeeping gene was calculated via BestKeeper V1. With the SD value less than 1, the housekeeping gene is considered as the stable one. Furthermore, the lower the SD value is, the higher stability that gene displays. Reversely, the gene was counted as being not stable if the SD value is higher than 1. Our results demonstrated that the U2AF gene exhibited the highest stability in total, particularly in the root sample treated with Hg and leaf samples treated with Pb, Hg, and As (Table 1). UCE was found to be the most stable reference gene in leaf samples treated with Cd, while eEF4a was identified to be the most stable reference gene for root samples treated with Cd and As. In the root and leaf samples treated with Pb and Cr, respectively, the UBC gene was shown to have the highest stability. When the root samples were treated with Cr, FTSH4 was found to be the most stable gene. Among all of the reference genes, the SD values of the CYP2 gene were demonstrated to be greater than 1 in all treatments except for PbR, indicating the instability of this gene (Table 1).

NormFinder Analysis
The lower value calculated from NormFinder V20 indicates the higher stability of the housekeeping gene expression. Results demonstrated that U2AF was shown to be the most stable reference gene in leaf samples under stress by Cd and root samples under stress by Hg, Cr, and As (Table 2). CYP5 was ranked as the most stable reference gene in root tissues treated with Cd. In response to Pb stress, 18S rRNA and eEF1a were identified to have the highest stability in both leaf and root tissues. In the leaf samples treated with Hg and Cr, the ACT12 gene was identified to have the highest stability. UBQ2 was shown to be the most stable reference gene in leaves treated with As. Overall evaluation identified four reference genes to be the most stable ones as follows: U2AF (0.463), CYP5 (0.484), UBQ1 (0.539), eEF4a (0.569). CYP2 was shown to be the least stable gene in all samples except for the root tissues treated with Pb and Hg ( Table 2).

RefFinder Analysis
RefFinder V1.0 was used to evaluate the comprehensive stability of reference genes integrating the methodologies of GeNorm, NormFinder, and BestKeeper analyses. Our results identified U2AF along with different genes to be the ideal reference genes in leaf tissues treated with Cd, Hg, Cr, As, and as well as root tissues in response to Pb, Hg, Cr, and As stress (Table 3). 18S rRNA and ACT12 were found to be optimal reference genes in Pb treated leaf tissues. In Cd treated root samples, CYP5 and UBQ1 were demonstrated to be the appropriate reference genes. Regarding the unstable reference genes, CYP2 was shown to have the least stability in all treatments except for root tissues with Pb. Furthermore, ACT2 was found to have poor stability in all samples except for root tissues treated with Pb and As (Table 3).

Detection of Target Gene Expression Levels Normalized by Screened Reference Genes
Recent studies have identified numerous genes involved in heavy metal stress response. The expressions of these genes were shown to function significantly in cellular transportation and enrichment of heavy metals as well as enhancement of plant resistance [18,19]. Especially gene coded proteins involved in the transportation of heavy metals have been extensively studied. These genes include the ABC transporter family, ZIP family, and heavy-metal ATPases (HMA) family, and they are reported to improve the tolerance towards heavy metals in plants [20][21][22]. Based on the comprehensive evaluation of reference genes stability in young leaves and young roots of switchgrass in response to different heavy metals, we selected and combined the reference genes with the least and highest stability to analyze the expression of PvZIP4 gene from ZIP family and PvPDR8 gene from ABC family in leaf and root tissues to validate our candidate reference genes (Figure 4). When CYP5/UBQ1 or U2AF/ACT12 was used as the reference gene pair, the change in relative expressions of target genes in response to Cd treatment was basically consistent. Meanwhile, the expression profile of target genes under Pb stress was also consistent when normalized by eEF1a/U2AF or 18S rRNA/ACT12, indicating that the expression of these reference genes is stable (Figure 4). However, when CYP2 was used as the reference gene, a significant decrease in the relative expression of PvZIP4 in Cd treated root samples shown at 3 h, 6 h, 12 h, 24 h, and 48 h compared to 0 h, was observed, which contradicted with the increasing trend demonstrated by other reference genes. Moreover, by using the CYP2 as the reference gene, the expression of PvPDR8 increased in leaves after 12, 24, and 48 h of Pb treatment compared to 1.5h, 3h, and 6h, as opposed to the time-dependent expression change pattern obtained from the reference genes 18S rRNA, ACT12, or the pair (18S rRNA/ACT12). Our results suggested that CYP2 had the poor stability of expression in young leaf tissues with different treatments, making it not appropriate to adjust the relative expression of target genes in an accurate manner.

Detection of Target Gene Expression Levels Normalized by Screened Reference Genes
Recent studies have identified numerous genes involved in heavy metal stress response. The expressions of these genes were shown to function significantly in cellular transportation and enrichment of heavy metals as well as enhancement of plant resistance [18,19]. Especially gene coded proteins involved in the transportation of heavy metals have been extensively studied. These genes include the ABC transporter family, ZIP family, and heavy-metal ATPases (HMA) family, and they are reported to improve the tolerance towards heavy metals in plants [20][21][22]. Based on the comprehensive evaluation of reference genes stability in young leaves and young roots of switchgrass in response to different heavy metals, we selected and combined the reference genes with the least and highest stability to analyze the expression of PvZIP4 gene from ZIP family and PvPDR8 gene from ABC family in leaf and root tissues to validate our candidate reference genes (Figure 4). When CYP5/UBQ1 or U2AF/ACT12 was used as the reference gene pair, the change in relative expressions of target genes in response to Cd treatment was basically consistent. Meanwhile, the expression profile of target genes under Pb stress was also consistent when normalized by eEF1a/U2AF or 18S rRNA/ACT12, indicating that the expression of these reference genes is stable (Figure 4). However, when CYP2 was used as the reference gene, a significant decrease in the relative expression of PvZIP4 in Cd treated root samples shown at 3 h, 6 h, 12 h, 24 h, and 48 h compared to 0 h, was observed, which contradicted with the increasing trend demonstrated by other reference genes. Moreover, by using the CYP2 as the reference gene, the expression of PvPDR8 increased in leaves after 12, 24, and 48 h of Pb treatment compared to 1.5h, 3h, and 6h, as opposed to the time-dependent expression change pattern obtained from the reference genes 18S rRNA, ACT12, or the pair (18S rRNA/ACT12). Our results suggested that CYP2 had the poor stability of expression in young leaf tissues with different treatments, making it not appropriate to adjust the relative expression of target genes in an accurate manner.

Discussion
The selection of proper internal reference genes largely depends on the target tissues, cells, and experimental conditions. The internal reference genes display specific differential expression when tissues/organs, developmental stages, and biotic/abiotic stress differ. There is no single internal reference gene which is constantly stable in expression under different experimental conditions in plants. The stability of the conserved reference genes also differs in different plant species. Take the internal reference gene GAPDH as an example, it displays poor expression stability in crops, such as wheat (Triticum aestivum L.) [23] and chicory (Cichorium intybus L.) [24], while being stable in grape (Vitis vinifera L.) [25] and sugarcane (Saccharum officinarum L.) [26]. Nakayama et al. [27] identified that Fbox/60s and Fbox/ABC are proper internal reference genes in the seedling tissue of soybean (Glycine max), while ELF1B and ACTB are identified as appropriate internal reference genes in soybean root tissues. It supports the argument that the internal reference genes display differential expression in different species and tissues. The previous studies showed that the expression of internal reference genes is closely related to experimental conditions. TIP41 in Arabidopsis thaliana L. was stable in expression under the nutrient deficient stress [28], while its expression stability was significantly reduced when the stress was induced by copper and cadmium [29]. In addition, the stability of the reference genes from the same gene family differs. Gutierrez et al. [30] discovered that the expression of the internal reference gene UBQ5 was more stable in Arabidopsis than UBQ4, UBQ10, and UBQ11 though they are from the same gene family. Therefore, it is imperative to select the stable reference gene based on the specific experimental conditions for RT-qPCR.
Switchgrass has gained increasing popularity in the study of energy plants in recent years. The previous studies focused on the expression analysis of internal reference genes in switchgrass under various abiotic stress conditions, such as drought, salinity, high temperature, and water flooding [10,12]. However, it did not delve into the systemic comparative analysis of expression stability of these internal reference genes among different tissues of switchgrass under heavy metal stress, resulting in the possible deviations in quantifications of the expression of target genes in response to different heavy metal treatments. Housekeeping genes are often expressed in cells with an active metabolism, which maintain the basic functions of cells and play an important role in the regulation of the cell cycle. Housekeeping genes are better candidate genes for evaluating gene expression and its functional properties [31]. Commonly used housekeeping genes in Gramineae crop, e.g., actin (ACT), 18S ribosomal RNA protein (18S rRNA), cyclophilin (CYP), eukariotic elongation factor (eEF), splicing factor U2af (U2AF), ftsh protease (FTSH), ubiqutin (UBQ), and ubiqutin conjugating enzyme (UCE) are considered to be suitable due to their presence in all nucleated cell types and essential functions in cell survival. Moreover, their expression has been considered to be stable in various tissues [32,33]. Our study reports the first validation of housekeeping genes in switchgrass allowing the identification of the most suitable reference gene(s) for normalization of RT-qPCR in different plant tissues (roots and leaves) and different time-courses subjected to heavy metal treatments such as by Cd, Pb, Hg, Cr, and As.
To avoid the limitations of using only single software analysis, our study applied three analytical approaches GeNorm, NormFinder, and BestKeeper to determine the expression stability of internal reference genes in different tissues under heavy metal stresses. The basis for evaluating gene stability in GeNorm is to use the logarithmic conversion value (2 −∆Ct ) of each gene to calculate the average variability (M value) [13]. Meanwhile, GeNorm can determine the optimal number of reference genes required for quantitative analyses: When the comparative difference analysis is performed on the internal reference gene normalization factor (V n /V n+1 ), the n value equals the number of optimal reference genes applied in RT-qPCR analysis. In this study, the gene expression analysis needs two reference genes to achieve the best performance. The combined use of the four house-keeping genes could be suitable for testing all the considered tissues, and stress conditions were based on the pairwise variation V n /V n+1 value. However, considering we commonly select only single heavy metal treatment to study the gene regulatory network of plants in response to heavy metal stress, four reference genes for gene normalization in all the considered tissues and stress conditions have limited practical values. The algorithm of NormFinder is similar to GeNorm using the logarithmic conversion value (2 −∆Ct ) as the relative expression of the gene to calculate the stability of gene expression [14]. BestKeeper focuses on the standard coefficient variation (SD) and variation correlation coefficient (CV) to screen the stability of internal reference genes [15]. The results from these approaches differ due to the algorithm differences. For example, U2AF was determined to be the most stable reference gene by GeNorm and NormFinder. However, it turned out that ACT12 was evaluated as the best in BestKeeper analysis. In addition, GeNorm alone suggested that eEF4a and U2AF were the top selections with the highest stability. The result inconsistency has been reported previously when different analytical software was applied [34]. Therefore, RefFinder integrates the algorithms of GeNorm, NormFinder, and BestKeeper to achieve comprehensive evaluation of the stability of reference genes, avoiding the unilateral judgment from a single method. In this study, U2AF displayed the highest stability in Cd, Hg, Cr, and As treated leaves as well as in Pb, Hg, Cr, and As treated roots of switchgrass. U2AF has been reported in other studies to be the most stable reference gene, such as in Pinus massoniana L., at different stages post-inoculation by nematode [35] and roots/leaves of Paspalum vaginatum Sw. under Cd and cold stress [36]. In the roots of switchgrass subjected to Cd treatment, CYP5 and UBQ1 were determined to be the most appropriate internal reference genes. CYP5 was identified as the most stable reference gene in ganoderma under various experimental conditions [37]. De Andrade et al. [38] found that UBQ1 was stably expressed in leaves of Saccharum spp under the drought stress. 18S rRNA displayed the best stability in Pb treated leaves of switchgrass.
Early selection of internal reference genes is mainly dependent on the assumptions of the essentiality of the housekeeping genes' functions. For example, based on the essential functions of Actin and TUB in cytoskeleton composition, these genes were speculated to be stably expressed in all cellular and physiological states [39]. However, the stability of internal reference genes can actually vary in different conditions in reality [40]. Actin2 in our study had poor stability under all experimental conditions. Some previous studies claimed that Actin was not the proper internal reference gene in chrysanthemum under abiotic stresses (high temperature, water flooding, aphid) [41] and wild type potato before and after the cold acclimation [42]. However, ACT12, another member of the Actin gene family, was shown to be stably expressed in leaves of switchgrass under the stress of Pb in our study. Thus, the expressions of housekeeping genes are not universally stable among various species in response to different stress conditions. The selection of internal reference genes should be expanded beyond the housekeeping gene families. Cyclophilins (CyPs) are ubiquitous proteins functioning in the folding of certain proteins involved in signal transduction processes. In Solanum tuberosum L., the level of a cyclophilin gene mRNA accumulation is stimulated by the application of abscisic acid and methyl jasmonate in tubers. However, treatment with fungal elicitor or salicilic acid has no such obvious effect [43]. Moreover, CYP mRNA synthesis was also shown to be variable in maize and bean with mercuric chloride treatment and in other abiotic stresses conditions, such as heat, wounding, salt stress, and low temperature [44,45]. In addition, different drugs significantly induced CYP transcription in human tissues [46]. In this study, CYP2 was found to be unstable as a reference gene in all samples except for root tissues treated with Pb. Therefore, particular caution should be taken when CYP is considered as a reliable reference gene.
To validate the selected internal reference genes from the screen, we chose PvZIP4 and PvPDR8 encoding the metal transporters as the target genes with expression induced by heavy metal stress [20,22]. The expression of these target genes in response to Cd and Pb stress was analyzed using two optimal and one poorly stable reference genes. Our results demonstrated that the target genes exhibited a general expression pattern in response to heavy metal stress when CYP5/UBQ1, U2AF/ACT12, eEF1a/U2AF, and 18S rRNA/ACT12 were used as the internal reference genes, while irregular patterns were shown with CYP2 selected to be the reference gene for RT-qPCR analysis. It suggested that the selection of proper internal reference genes is essential for RT-qPCR quantitative analysis.

Conclusions
Our study provides a good reference for selecting proper internal reference genes to study the expression of target genes under heavy metal stress in switchgrass. According to the results of RefFinder, U2AF was the best reference gene in Cd, Hg, Cr, and As treated leaves as well as in Hg, Pb, As, and Cr stressed root tissues. 18S rRNA was considered the most stable reference gene in Pb treated leaf tissues. CYP5 was determined to be the optimal reference gene in Cd treated root tissues. While the least stable reference gene was identified to be CYP2 in all tested samples except for root tissues stressed by Pb. In addition, the choice of the combination of the appropriate internal reference genes can significantly impact the analysis of the target gene expression pattern in response to different heavy metal stresses. Our findings established a solid foundation to further study the gene regulatory network of switchgrass in response to heavy metal stress.

Conflicts of Interest:
The authors declare no conflict of interest.