Selection and Validation of Reference Genes for qRT-PCR Gene Expression Analysis in Kengyilia melanthera

Kengyilia is a newly established genus. Most species in this genus survive in hash environment, which might be an indicator of an acquirement of stress resistance genes and the potential for molecular breeding in Triticeae species. Quantitative real-time PCR (qRT-PCR) is a widely used technique with varied sensitivity heavily dependent on the optimal level of the reference genes. K. melanthera is a typical psammophyte species which has high drought resistance. The reference genes of K. melanthera are not yet reported. This study aims to evaluate the expression stability of 14 candidate reference genes (EF1A, GAPDH, ACT1, UBI, TUBB3, TIPRL, CACS, PPP2R1B, TUBA1A, EIF4A1, CYPA3, TCTP, ABCG11L, and FBXO6L) under five treatments (drought, heat, cold, salt, and ABA) and find the most stable and suitable one even upon stressed conditions. The software NormFinder, GeNorm, BestKeeper, and RefFinder were used for data analysis. In general, the genes CACS and PPP2R1B are concluded to have the best overall performance under the various treatments. With the ABA treatment, TCTP and TIPRL show the best stability. CACS and TCTP, as well as TIPRL and CYPA3, were most stable under the treatments of cold and salt, respectively. CACS and FBXO6L were ranked the highest with the heat treatment and drought treatment, respectively. Finally, the Catalase-1 (CAT1) gene was used to verify the reliability of the above reference genes. Accordingly, CAT1’s expression pattern remained unchanged after normalization with stable reference genes. This study provides beneficial information about the stability and reliability of potential reference genes for qRT-PCR in K. melanthera.


Introduction
There are about 26 species and six varieties of Kengyilia, a genus in the tribe Triticeae established in 1990 [1]. K. melanthera is a self-pollinated perennial grass species with allohexaploid. It is a typical psammophyte species, mainly distributed in sandy river banks, dunes, and sandy meadows of the Qinghai-Tibet Plateau at an altitude of 3300-4750 m [2]. As such, K. melanthera has strong drought and cold resistance. It can be used as a grass for plateau desertification control and ecological restoration, or as plateau forage because of its high biomass and nutritional quality. In addition, some studies reported that K. melanthera has a strong resistance to wheat head scab [3]. Kengyilia and other Triticeae perennials have a vast genetic reservoir, which might be used to improve annual cereals [4]. Wheat (Triticum aestivum) is one of the major staple crops in the world, of which the yield is limited by abiotic stress and biotic stress [5]. Among cultivated wheat, the breeding potential is already exhausted as there is an increasingly narrow range of genetic variation [6]. At present, certain genes of related species are transferred to wheat to improve its yield and stress resistance [5,7,8]. We consider it necessary to explore K. melanthera's resistance genes and related regulatory pathways, given their potential for wheat improvement.
Gene expression analysis is a useful tool in exploring the mechanisms of stress resistance in K. melanthera. Quantitative real-time polymerase chain reaction (qRT-PCR) is a frequently used technique to examine the gene expression level with the strength in fast reaction, precision, specificity, sensitivity, and repeatability [9][10][11][12]. Absolute and relative quantification are two methods of quantitative gene expression [13,14]. However, it is usually unnecessary to express quantitative data as an absolute copy number [14]. Therefore, relative quantification is more commonly applied [13]. In most cases, the reliability of qRT-PCR is affected by multiple factors, including RNA quantity and quality, the efficiency of cDNA synthesis, the quantity of the starting DNA template, and so on [10,12,15]. Hence, to compensate for the above difference, it is inevitable to have a requirement for reliable reference genes for normalization.
The ideal reference genes are supposed to exhibit the least expression level difference under different environmental stresses or in different tissues, organs, and developmental stages [9,11,15]. In qRT-PCR analysis for plant species, common reference genes include Actin, Polyubiquitin (UBQ), Elongation factor-1a (EF-1α), α-Tubulin (TUA), 18S rRNA, β-Tubulin (TUB), Cyclophilin (CYP), and Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) [16][17][18]. However, there are quite a few studies reporting that the expression level of these reference genes is not stable under different conditions [19][20][21][22]. Inappropriate selection of the unsuitable reference genes can lead to significant mistakes in results [23]. Therefore, identifying ideal reference genes for a particular species and understanding their performance under different conditions are critical for delivering reliable gene expressional analytical results [13,24].
Currently, no reference genes have been described for K. melanthera. Herein, the expression stability of 14 reference gene candidates in leaves under different experimental situations (drought, heat, cold, salt, and abscisic acid) was identified and verified by comparison to the Catalase-1 (CAT1) gene. This work aims to contribute to the future exploration and utilization of K. melanthera's resistance genes.

Plant Materials and Treatments
K. melanthera seeds were kindly provided by Sichuan Academy of Grassland Sciences (Chengdu, China). Seeds were planted in plastic pots (20 × 15 × 5 cm) with quartz sand. A total of 1.5 g of seeds were used for each pot. Hoagland nutrient solution was applied. The pots were placed in a growth chamber at 25 • C/20 • C (day/night). A 12 h photoperiod was set for the growth chamber. Twenty-one-day old plants were used for all experiments. Five treatments were set up for the screening of the reference genes. Each treatment has three biological replicates. In abscisic acid (ABA) treatment, the Hoagland nutrient solution was added with 50 µM ABA. In cold and heat treatment, the plants were placed in an incubator with a temperature setting of 4 • C or 38/33 • C (day/night), respectively. In salinity treatment, the Hoagland nutrient solution supplemented with 250 mM NaCl was applied. In drought treatment, PEG6000 at 20% concentration was performed. Finally, the leaf samples were collected at five time points, respectively-0, 12, 24, 48, and 96 h post treatment. The harvested samples were stored in a minus 80 • C lab freezer.

RNA Extraction and cDNA Synthesis
The RNA extraction kit (Monad Biotech, Suzhou, China) was used for the total RNA extraction following the manufacturer's instructions. BioPhotometer (Eppendorf, Hamburg, Germany) measurement and 1% agarose gel were conducted to examine the RNA concentration and integrity. An absorbance ratio of A260/A280 ranging from 1.8 to 2.2 and A260/A230 ratio equal to 2.0 were preferred. Based on the recommendation of the Evo M-MLV RT Mix Kit (Accurate Biotech, Changsha, China), 0.8 µg total RNA was measured to perform the cDNA synthesis.

Primer Design
Fourteen reference gene candidates and one target gene were selected and named according to the sequence similarities to known genes (Table 1). These genes were obtained from our full-length transcriptome data (accession number: PRJNA735213) of K. melanthera by BLAST search using reported gene sequences. The Primer-BLAST tool by NCBI (https: //www.ncbi.nlm.nih.gov/tools/primer-blast/, (accessed on 1 May 2021)) was used for primer design. The principles of primer design were as follows: annealing temperature at 58-62 • C (optimal Tm was 60 • C), primer length at 18-25 bp, GC content at 40-60%, and the length of amplification product between 80-200 bp. Primer specificity was detected by conventional PCR and 2% agarose gel electrophoresis.

Quantitative RT-PCR Amplification
CFX96 PCR detection system (Bio-RAD, Hercules, CA, USA) was used to perform qRT-PCR. The 10 µL final reaction volume included 5 µL 2 × SYBR Green Premix (Monad Biotech, Suzhou, China), 1 µL cDNA, 1 µL of forward and reverse primer (final concentration of 0.2 µmol·L −1 ), and 3 µL ddH 2 O. The PCR program was set as following: pre-denaturation at 94 • C for 20 s, followed by 40 cycles of denaturation at 94 • C for 10 s, and annealing at 60 • C for 20 s. To verify primer specificity, T m and melting curves were analyzed between 65 • C-95 • C, and fluorescence values were detected at each 0.5 • C increase. All samples used for qRT-PCR analysis were set with three technical replicates. Standard curves for each candidate reference gene were constructed to determine PCR amplification efficiency and regression coefficients (R 2 ), and the data were further analyzed by CFX Manager Software 3.1.

Analysis of Reference Gene Candidates' Expression
The stability of 14 candidate reference genes was evaluated by GeNorm [25], NormFinder [26], BestKeeper [27], and online tool RefFinder [28]. Prior to GeNorm and NormFinder analysis, the raw Cq values need to be converted into relative quantities according to the formula 2 −∆Cq (∆Cq = each corresponding Cq value − lowest Cq value) [29]. GeNorm determined the stability of reference genes by calculating M values; the reference genes with better stability have smaller M values. GeNorm can also calculate pairwise variation (V). When the V n /V n+1 value ≤ 0.15, the number of suitable reference genes is "n" [25]. NormFinder selected the most suitable reference gene by calculating the stability value of candidate gene expression. A low stability value indicates that the gene is stable [26]. BestKeeper calculated standard deviation (SD) and coefficient of variation (CV) by raw Cq values. The more stable reference gene had the lower SD value [27]. RefFinder can generate a comprehensive ranking based on the analysis results of GeNorm, BestKeeper, NormFinder, and Delta Ct.

Validation of Reference Genes
The two most stable reference genes, alone and in combination, and the least stable reference gene were used to normalize the expression of the target gene, namely CAT1. The results were calculated by 2 −∆∆Cq method [30]. SPSS 22 (IBM, Armonk, NY, USA) was used for statistical significance analysis.

Primer Specificity and Amplification Efficiency
In this study, the sequences of 14 reference gene candidates and one target gene were extracted from the full-length transcriptome of K. melanthera. All had >90% identity with the sequences of homologous genes registered in the NCBI, and their E values were zero (Table S1). The results of PCR amplification and 1% agarose gel electrophoresis (matching the expected amplicon size and single amplicon fragments on the gel) further proved the sequence accuracy ( Figure S1). In addition, the melting curve of each pair of primers had only one peak, indicating their specificity ( Figure S2). The primer amplification efficiency was in the range of 91.1% and 106.5%, and the regression coefficients (R 2 ) was in the range of 0.985 and 0.999 (Table 1). Therefore, our qRT-PCR data was reliable.

Expression Profile of the 14 Reference Gene Candidates in Response to Different Treatments
The expression abundance of 14 reference gene candidates in all samples are demonstrated in Figure 1 These results suggested that these candidate genes were not stably expressed in the stressed conditions. It is imperative to conduct a screening to find out appropriate reference genes in K. melanthera under certain conditions.

Analysis of Reference Genes Stability
A total of 14 reference candidate genes are subjected to the evaluation of the expression stability in response to five treatments with the help of NormFinder, BestKeeper, and GeNorm. The final overall ranking was calculated by RefFinder. GeNorm also calculated the paired variation (V n/n+1 ) value to determine the optimal number of reference genes for the normalization of target gene expression. If the V n/n+1 value is <0.15, the optimal number of reference genes is n; if not, another reference gene should be brought into consideration [25]. In this study, V 2/3 values were <0.15 under each treatment alone, while V 2/3 and V 3/4 values were >0.15 and V 4/5 values were <0.15 in all samples ( Figure 3). The results showed that it was adequate to use the combination of two reference genes in a separate treatment, but the combination of four reference genes was needed in the synthesis of all treatment samples.

NormFinder Analysis
The expression stability values of 14 reference gene candidates were calculated and ranked by NormFinder (Figure 4). The lower the stability value indicates the higher stability of the reference gene expression. In cold, heat, drought, and all samples, CACS was ranked as the most stable reference gene, with stability values of 0.005, 0.165, 0.164, and 0.182, respectively. Under ABA treatment, TCTP was the most stable gene, with a stability value of 0.141. TIPRL was the most stable gene under salt treatment, with a stability value of 0.142. In various treatments, TUBA1A and ABCG11L were ranked lowest, indicating their poor stability, consistent with the results of GeNorm analysis.
nes 2022, 13, x stressed conditions. It is imperative to conduct a screening to find out ap ence genes in K. melanthera under certain conditions.

Analysis of Reference Genes Stability
A total of 14 reference candidate genes are subjected to the evaluatio sion stability in response to five treatments with the help of NormFinder, GeNorm also calculated the paired variation (Vn/n+1) value to determine the optimal number of reference genes for the normalization of target gene expression. If the Vn/n+1 value is <0.15, the optimal number of reference genes is n; if not, another reference gene should be brought into consideration [25]. In this study, V2/3 values were <0.15 under each treatment alone, while V2/3 and V3/4 values were >0.15 and V4/5 values were <0.15 in all samples (Figure 3). The results showed that it was adequate to use the combination of two reference genes in a separate treatment, but the combination of four reference genes was needed in the synthesis of all treatment samples.   GeNorm also calculated the paired variation (Vn/n+1) value to determin number of reference genes for the normalization of target gene expressio value is <0.15, the optimal number of reference genes is n; if not, another r should be brought into consideration [25]. In this study, V2/3 values were <0. treatment alone, while V2/3 and V3/4 values were >0.15 and V4/5 values we samples ( Figure 3). The results showed that it was adequate to use the comb reference genes in a separate treatment, but the combination of four referen needed in the synthesis of all treatment samples.   bility of the reference gene expression. In cold, heat, drought, and all samples, CACS was ranked as the most stable reference gene, with stability values of 0.005, 0.165, 0.164, and 0.182, respectively. Under ABA treatment, TCTP was the most stable gene, with a stability value of 0.141. TIPRL was the most stable gene under salt treatment, with a stability value of 0.142. In various treatments, TUBA1A and ABCG11L were ranked lowest, indicating their poor stability, consistent with the results of GeNorm analysis.

BestKeeper Analysis
BestKeeper calculated the CV and SD according to the original Cq value [27]. The lower the CV ± SD, the more stable the gene expression. The most stable genes under ABA, cold, and drought treatment were EF1A, TIPRL and PPP2R1B, respectively. The FBXO6L gene was identified to have the highest stability in heat and salt treatment, while its stability in ABA, cold stress, and all samples ranked in the middle. Obviously, the stability of a gene may change under different stress treatments. Among all samples, TIPRL was the highest stable gene. The stability of the ABCG11L gene ranked lowest of all treatments ( Table 2).

BestKeeper Analysis
BestKeeper calculated the CV and SD according to the original Cq value [27]. The lower the CV ± SD, the more stable the gene expression. The most stable genes under ABA, cold, and drought treatment were EF1A, TIPRL and PPP2R1B, respectively. The FBXO6L gene was identified to have the highest stability in heat and salt treatment, while its stability in ABA, cold stress, and all samples ranked in the middle. Obviously, the stability of a gene may change under different stress treatments. Among all samples, TIPRL was the highest stable gene. The stability of the ABCG11L gene ranked lowest of all treatments ( Table 2).

Validation of The Reference Genes Identified from This Study
CAT helps protect cells from H 2 O 2 by breaking it down into dioxygen and water. The expression levels of the CAT gene in plants is influenced by multiple [31]. To identify the reliability of the stability ranking of reference gene candidates, the relative expression levels of the CAT1 gene were normalized using the two highest and least stable genes in the RefFinder ranking. Except for 48 h of drought treatment, there is no significant difference in CAT1 expression using the two most stable genes as reference genes at all other time points (p < 0.05) ( Figure 5). When the least stable gene was used as reference, CAT1 expression was significantly different in most cases (p > 0.05). For each treatment, the expression pattern of CAT1, analyzed by the two stable reference genes, remained unchanged. Under cold treatment, CAT1 expression showed an increasing trend ( Figure 5B); under heat, it showed a decrease-increase-decrease change ( Figure 5C); under ABA and salt and drought treatment, CAT1 first showed an increase and then a decrease ( Figure 5A,D,E). However, when normalized by the most unstable reference genes, CAT1 expression patterns tended to be abnormal. These results showed that using unstable genes used the expression results of target genes as reference biases.

Discussion
Kengyilia is a newly established plant genus [1]. Most species in this genus grow in harsh environments [2]; therefore, having excellent stress resistance genes of great potential value for Triticeae plant breeding. qRT-PCR plays a key role in gene expression research. Appropriate reference genes greatly enhance the precision and consistency of qRT-PCR in plants [32]. In the past, most reports focused on Kengyilia's classification and phylogeny [1,[33][34][35], but there was no report on the screening of reference genes. Currently, only a few plant genomes are published, excluding K. melanthera. Furthermore, compared to the second-generation transcriptome sequencing, the third-generation full-length transcriptome sequencing holds more advantages, including longer read length and less sequencing and assembly error, which is more beneficial for the study of plants without a reference genome [36,37]. In this study, 14 reference gene candidates were used and in-

Discussion
Kengyilia is a newly established plant genus [1]. Most species in this genus grow in harsh environments [2]; therefore, having excellent stress resistance genes of great potential value for Triticeae plant breeding. qRT-PCR plays a key role in gene expression research. Appropriate reference genes greatly enhance the precision and consistency of qRT-PCR in plants [32]. In the past, most reports focused on Kengyilia's classification and phylogeny [1,[33][34][35], but there was no report on the screening of reference genes. Currently, only a few plant genomes are published, excluding K. melanthera. Furthermore, compared to the second-generation transcriptome sequencing, the third-generation fulllength transcriptome sequencing holds more advantages, including longer read length and less sequencing and assembly error, which is more beneficial for the study of plants without a reference genome [36,37]. In this study, 14 reference gene candidates were used and investigated in order to find the most stable and appropriate ones in response to different treatments.
Here, three types of software (GeNorm, NormFinder, and BestKeeper) were used to evaluate the stability of 14 reference gene candidates under five different treatments. The stability of some genes was similar in the three platforms: for instance, TCTP had the highest while TUBA1A and ABCG11L had the lowest stability under ABA treatment; ABCG11L, UBI, and TUBA1A had the worst stability under cold treatment; and TUBA1A, ABCG11L, and TUBB3 had the worst stability under heat treatment. Nevertheless, the stability of some genes differed. Under drought treatment, the highest stable genes were FBXO6L/TCTP, CACS, and PPP2R1B, respectively. This may be owning to the algorithm difference among these three software [10].
According to RefFinder's results, the highest stable reference genes under various experimental situations may differ. TCTP was the highest stable gene under ABA treatment but ranked fifth under salt treatment. TIPRL was the highest stable gene under salt treatment; however, it ranked seventh in drought treatment and tenth in cold treatment. CACS was the highest stable gene under cold, heat, and drought treatment, but it ranked sixth in salt treatment. This is in agreement with other studies: for instance, Lycoris aurea's highest stable genes under drought and cold treatment were PTBP1 and CYP2, respectively, while EXP1 was the most stable gene under salt, heat, and ABA treatment [38]. Hence, it is crucial to conduct screening for the appropriate reference genes under different conditions for each species [39,40].
It Is worth noting that, although CACS is not a commonly used reference gene, it was the most stable under cold, heat, and drought treatment, of all samples in the RefFinder ranking of this study. CACS are subunits of the clathrin adaptor protein (AP) complex, through which the AP complex can collect cargo and recruit other proteins involved in the formation of clathrin-coated vesicles [41]. The CACS gene was reported as the highest stable reference gene in watermelons (Citrullus lanatus) under abiotic stress (cold, salt, drought) [42], and in cucumber (Cucumis sativus) under both short-term and long-term heavy metal stress (Cd, Pb, Cu, Zn, Mn, Ni) [43]. TUB is a traditional reference gene, but its stability ranked last in this study [18]. Abnormalities also appeared when TUBB3 was used as a reference gene to test the expression pattern of CAT1 gene. These data suggests that conventionally and commonly used reference genes may not be a good choice for all species and all experimental conditions. New reference gene candidates can be further explored.
Some studies have shown that normalization of target gene expression using one reference gene may lead to large experimental errors [25], which can be reduced by using ≥2 reference genes. According to GeNorm analysis, V 2/3 values were <0.15 in all treatments; therefore, only two reference genes were enough under a single treatment condition. When considering all samples, the V 2/3 and V 3/4 values were >0.15, and the V 4/5 value was <0.15; hence, four reference genes were recommended. However, in the final validation, we found that stable reference genes, either alone or in combination, could also correctly reveal the expression pattern of the CAT1 gene. Therefore, although ≥2 reference genes can increase the precision of the expression analysis, the number of reference genes depends on experimental conditions [15,40].

Conclusions
This study is the first report that reveals a series of appropriate reference genes for K. melanthera under different experimental conditions. The highest stable reference genes were CACS and PPP2R1B in all samples; TCTP and TIPRL under ABA treatment; CACS and TCTP under cold treatment; CACS and FBXO6L under both heat and drought treatment; and TIPRL and CYPA3 under salt treatment, respectively. This study greatly facilitates the gene expression analysis of K. melanthera under different experimental conditions. It also paves the way for the screening for reference gene candidates in other related species.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13081445/s1, Figure S1: Amplification products of the fourteen candidate reference genes and target gene.; Figure S2: Melting curves of the fourteen candidate reference genes and target gene; Table S1: Description of 14 candidate reference genes and one target gene.