Reference Gene Selection for Quantitative Real-Time PCR of Mycelia from Lentinula edodes under High-Temperature Stress

Housekeeping genes are important for measuring the transcription expression of functional genes; 10 traditional reference genes, TUB, TUA, GADPH, EF1, 18S, GTP, ACT, UBI, UBC, and H2A, were tested for their adequacy in Lentinula edodes (L. edodes). Using specific primers, mRNA levels of these candidate housekeeping genes were evaluated in mycelia of L. edodes, which were treated with high-temperature stress at 37°C for 0, 4, 8, 12, 18, and 24 hours. After treatment, expression stability of candidate genes was evaluated using three statistical software programs: geNorm, NormFinder, and BestKeeper. According to geNorm, TUB had the lowest M values in L. edodes strains 18 and 18N44. Using NormFinder, the best candidate reference gene in strain 18 was TUB (0.030), and the best candidate reference gene in strain 18N44 was UBI (0.047). In BestKeeper analysis, the standard deviation (SD) values of UBC, TUA, H2A, EF1, ACT, 18S, and GTP in strain 18 and those of GADPH and GTP in strain 18N44 were greater than 1; thus, these genes were disqualified as reference genes. Taken together, only UBI and TUB were found to be desirable reference genes by BestKeeper software. Based on the results of three software analyses, TUB was the most stable gene under all conditions and was verified as an appropriate reference gene for quantitative real-time polymerase chain reaction in L. edodes mycelia under high-temperature stress.


Introduction
Quantitative real-time polymerase chain reaction (qRT-PCR), which is characterized by high sensitivity, strong specificity, a repetitive dynamic quantitative range and high throughput, is one of the most commonly used techniques for gene expression analysis. The accuracy of qRT-PCR results is largely dependent on the selected reference genes [1][2][3][4][5][6], the validity of which is a prerequisite for the correct application of qRT-PCR to analyze changes in target gene expression [7][8][9]. To obtain more accurate and reliable results, reference genes are required for standardized measurement of the expression levels of target genes [10]. In recent years, considerable research has shown that the choice of reference genes should be associated with the experimental conditions of the investigator. The selection of reference genes may change depending on experimental conditions [1,2,[11][12][13]. In general, a basic component of the cytoskeleton or a relatively stable expressed gene involved in basic metabolic processes is selected as the reference gene. Nonetheless, research has shown that the expression levels of such stably expressed genes may not be stable in various cell types, at different developmental stages or under different experimental conditions [14][15][16][17]. Therefore, the evaluation and selection of reference genes under different experimental conditions are of great significance for obtaining reliable and precise qRT-and edible fungus [21,22]. This mushroom fructifies under moderate and low temperatures, with an optimal temperature for spore germination between 22 and 26 ∘ C, for mycelial growth between 24 and 27 ∘ C, for primordium differentiation between 10 and 12 ∘ C, and for fruit body development between 8 and 20 ∘ C [22]. If the temperature is below its optimal value for mycelial growth, L. edodes will grow more slowly, and mycelia will be thicker. Poor growth may occur both below 0 ∘ C and above 32 ∘ C, whereas growth cessation and terminal death occur at 35 ∘ C and 38 ∘ C, respectively. If the temperature is above the optimum, the fruit body grows faster, with a thin cap and a long stipe that result in a loose texture, an easyto-open umbrella, and poor quality. In this study, L. edodes mycelia were subjected to high-temperature stress for various time periods. Ten genes were selected as candidate reference genes, and three types of software (geNorm, NormFinder, and BestKeeper) were used to evaluate these candidates as reference genes for L. edodes under high-temperature stress. The results of this study lay the foundation for future research regarding molecular mechanisms of L. edodes in response to high-temperature stress.

Material and Processing. L. edodes strains 18 and 18N44
were obtained from the Institute of Edible Fungi, Shanghai Academy of Agricultural Sciences.
Mycelia of L. edodes were incubated and then treated with high-temperature stress at 37 ∘ C for 4, 8, 12, 18, and 24 hours; individuals not exposed to high-temperature stress were used as the control group. Biology duplications are three times. The mycelia were collected under sterile conditions, flash-frozen in liquid nitrogen, and stored at -80 ∘ C.

Extraction and Purification of Mycelia Total RNA.
Total RNA of mycelia was extracted with a Redzol reagent kit from Beijing SBS Genetech Co., Ltd., according to the instruction manual. RNA integrity was assessed by 1% agarose gel electrophoresis at 120 V for 10 minutes, and RNA purity was measured using a NanoDrop 2000 C (Thermo Fisher Scientific, Hudson, New Hampshire, USA).
The primers were designed using the NCBI online server tool with an amplification product length ranging from 70 to 300 bp. All primers were synthesized by the Shanghai Generay Biotech Co., Ltd., and the primer sequences are presented in Table 1.

Evaluation of Amplification Efficiency of Candidate Reference
Genes. Using L. edodes mycelia that were not treated with high-temperature stress (hereafter referred to as CK) for comparison, the primer pairs were amplified to detect amplicon production using cDNA as the template. The successfully amplified bands products were retrieved and ligated to a Tvector, and genes were sequenced by Sangon Biotech.
The amplification efficiency of the reference genes was calculated referring to the method of Li Peng [27]. The cDNA template of CK was diluted, and the cycle threshold (Ct value) was determined based on qRT-PCR. Calibration curves were prepared based on logarithmic fitting with the cDNA dilution and Ct value as the abscissa and ordinate, respectively. The amplification efficiency of each primer pair was calculated according to relevant data, with three repeats for each specimen.

Reverse Transcription (RT) of cDNA and Quantitative Real-Time (qRT) PCR.
This experiment was performed on ice using Prime Script RT Reagent Kit with gDNA Eraser (Takara, Shiga, Japan) according to the manufacturer's instructions. cDNA acquired by reverse transcription was stored at -20 ∘ C. qRT-PCR was performed using the SYBR5 Premix Ex Taq6 II (Takara, Shiga, Japan) Kit and a StepOnePlus Real-Time PCR instrument. The reaction system included 10 L SYBR5 Premix Ex Taq6 II (Tli RNaseH Plus) (2x), 0.4 L ROX dye, 2 L template cDNA, 0.4 L each primer, and 6.8 L RNase-free water. The PCR amplifications were performed using a 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) using the following: predegeneration at 95 ∘ C for 30 s; PCR at 95 ∘ C for 5 s, 60 ∘ C for 15 s, and 72 ∘ C for 15 s for 40 cycles; melting curve at 95 ∘ C for 15 min, 60 ∘ C for 30 s, and 95 ∘ C for 15 min.
qRT-PCR was performed on cDNA from hightemperature-sensitive and high-temperature-tolerant stains 18 and 18N44 under heat stress treatment and Ct values were recorded. The data were later imported into geNorm software, and the minimum value of Ct under different treatments was set to 1. Hence, the relative expression of other reference genes should be as follows: 2 -ΔΔCt (ΔCt = Ct-1, Ct value of each reference gene) [29]. The threshold M was set to less than 1.5, and NormFinder and BestKeeper were used for additional analyses.

Validation of Reference Genes.
In order to confirm the validity of the selected reference genes for data normalization, the expression profiles of 10 candidate reference genes were synthetically analyzed and sorted. Two most stable and one least stable candidate reference genes were selected. In the meanwhile, expression levels of hydrophobin gene (hyd1) and heat shock protein 100 gene (hsp100) under high-temperature stress were verified by qRT-PCR analysis.
The qRT-PCR amplification conditions were the same as those described above. The relative expression levels of the two genes were calculated by 2 -ΔΔCt method [30,31]. The analysis of gene expression by qRT-PCR was performed using three technical replicates for each biological replicate.

Total RNA Quality.
In this study, total RNA of L. edodes mycelia was extracted using the Redzol reagent method and  Figure S1, the 28S, 18S, and 5S subunit bands were clear and the fragments were intact, with no obvious degradation; therefore, the material was used in follow-up experiments.

Amplification Efficiency and Specificity of Candidate
Reference Gene Primers. PCR was performed to confirm the amplification capacity of all candidate reference gene primers. The results indicated a single band for all 10 pairs of primers; the amplified products were in line with the size of the expected fragment, and no primer dimers appeared. The PCR products were purified, and the sequencing results were consistent with reference sequences, which also found no similarity with other genes by Nucleotide-BLAST in the whole-genome shotgun contigs (WGS) database for strain W1-26. Verified by gradient PCR, all 10 pairs of primers showed excellent amplification at an annealing temperature approximately 60 ∘ C.

Calculation of Amplification Efficiency of Candidate Reference Genes.
The cDNA from CK samples were diluted into five concentration gradients, with each sample repeated three times. qRT-PCR was performed for all 10 candidate reference genes. There were no impurity peaks in the dissociation curve, i.e., only the main peak was observed (Supplementary Figure S2). The Ct values of the candidate reference genes were evaluated using the dilution concentration of cDNA as the abscissa and the Ct value as the ordinate on the logarithmic fitting calibration curve. Thus, the regression coefficients of reference genes were calculated in accordance with the formula = 10 -1/slope -1, and amplification efficiency can be determined based on the slope of the curve. The results showed that the amplification efficiency of each candidate reference gene fell between 95% and 110%. The value of R 2 was greater than 0.99, indicating a single dissolution curve. qRT-PCR was therefore performed using the ΔΔCT method ( Table 2).

geNorm Analysis Results.
For geNorm analysis, we followed the methods of Estchmann and others [32]. The minimum value of Ct of a specific reference gene was set to 1, and the relative expression of other reference genes was expressed using the 2 -ΔCt formula. By calculating the Ct values of the 10 reference genes of samples, M can be determined using geNorm. When the threshold value of M is smaller than 1.5, the candidate reference parameters are applicable. A stable expression map of candidate reference genes from two strains of L. edodes under high-temperature stress over different time periods is illustrated in Figure 1. The M value of the 10 genes in 18 and 18N44, including TUB, TUA, GADPH, EF1, 18S, GTP, ACT, UBI, UBC, and H2A, was determined using geNorm software. Based on the results, UBC at 2.742 in strain 18 and GTP at 1.508 in strain 18N44 exceeded the threshold, whereas the other reference genes were comparatively stable. The stability of these genes was evaluated using material from strains 18 and 18N44. With regard to strain 18, ACT and TUB showed high stability, followed by UBI and 18S; in contrast, UBC exhibited low stability (Figure 2(a)). Differences in the stability of the reference genes in strain 18N44 were also examined, with TUB, TUA, and UBI showing the most stable expression (Figure 2(b)). TUB displayed the lowest M values in both strains of L. edodes.  When the threshold value of M is smaller than 1.5, the candidate reference parameters are applicable. The two less stable genes (UBC and GTP) are eliminated.

NormFinder Analysis
Results. geNorm software can not significantly differentiate reference genes that share similar expression patterns. However, NormFinder software can effectively avoid unfavorable analysis results in this situation. NormFinder software conducts expression stability analysis on candidate reference genes mainly based on the results of variance analysis, and the software then lists genes in accordance with their stability [33]. The analysis of the expression stability of the 10 candidate reference genes in L. edodes at 37 ∘ C for different durations is presented in Table 3; lower values indicate more stable expression. The results showed that the best candidate reference gene in strain 18 was TUB (0.030); that in strain 18N44 was UBI (0.047).

BestKeeper Analysis
Results. The standard deviation (SD) and coefficient of variation (CV) as well as a pair of correlation coefficients (Poisson correlation coefficient) were calculated by BestKeeper via paired correlation analysis. Reference genes with an SD below 1 are considered stably expressed, and a smaller SD indicates a more stable reference gene [34]. BestKeeper analysis (Figure 3) revealed SD values greater than 1 for UBC, TUA, H2A, EF1, ACT, 18S, and GTP in strain 18 and GADPH and GTP in strain 18N44, disqualifying them as reference genes. Hence, only UBI and TUB were found to be desirable reference genes.

Determination of the Optimal Number of Reference Genes
for Normalization. Some reports have suggested that the use of more reference genes may lead to more stable results [13,35]. To determine the optimal number of reference genes in each experimental condition, pairwise variation (V) was calculated using geNorm by applying a cut-off value of 0.150 [13]. The pairing variable value Vn/n+1 can be analyzed using geNorm software to determine the optimal reference gene under certain conditions. In this study, Vn/n+1 was V3/4<0.150 (0.136, Figure 4(a)) in strain 18, indicating that three reference genes are sufficient as a standardized indicator. The variable value satisfied the condition of V2/3<0.150 (0.097, Figure 4(b)) in strain 18N44, suggesting that two reference genes should be used for normalization. According to NormFinder software analysis, the most stable reference gene in strain 18N44 is UBI, followed by TUA and TUB; in strain 18, the most stable reference gene is TUB, followed by EF1 and UBI. Analysis using BestKeeper (Table 4) software indicated that smaller SD and CV values were correlated with better outcomes. If the SD value exceeded 1, it was insufficient as a reference gene; hence, the eligible candidate reference genes were GADPH, TUB, and TUA. Based on the above analyses, TUB was chosen as the most suitable reference gene.

Validation of Reference Genes.
To validate the selection of candidate reference genes, the relative expression of two key genes hsp100 and hyd1 in response to high-temperature stress were monitored. The results revealed that the two genes were differentially or specifically expressed ( Figure 5).
The above comprehensive analysis indicates that, in the strain 18, TUB and ACT were the most stable candidate reference genes, and UBC was the least stable candidate reference gene. From Figures 5(a) and 5(b), relative expression levels of hsp100 and hyd1 obtained by normalizing against TUB or ACT indicated similar trends. However, normalization against UBC led to the opposite conclusion, exhibiting that expression of the two genes was relatively higher under hightemperature stress at 37 ∘ C for 8, 18, and 24 hours.  In the strain 18N44, TUB and TUA were the most stable candidate reference genes, and GTP was the least stable candidate reference gene. From Figures 5(c) and 5(d), relative expression levels of hsp100 and hyd1 obtained by normalizing against TUB or TUA indicated similar trends. However, normalization against GTP led to the opposite conclusion, exhibiting that the two genes were relatively higher under high-temperature stress at 37 ∘ C for 8, 18, and 24 hours.

Discussion
In qRT-PCR analysis, certain differences in the expression stability of the same reference gene can be observed among various species. Therefore, prior to expression analysis of a target gene, it is necessary to evaluate the expression stability of reference genes under certain experimental conditions [11,33,36]. The ideal reference genes should have unaltered expression levels during the entire experimental process and with proper expression intensity [3]. However, as there is no single reference gene that is appropriate for all experimental conditions [37], it is necessary to examine the expression stability of reference genes at different stages [38][39][40] or under different biotic and abiotic stresses to assess their use for normalization [41][42][43].
A stable reference gene is a prerequisite for improving the reliability of qRT-PCR results. The selection of reference genes was once mainly based on the functions of housekeeping genes. For instance, Actin and TUB are the basic components of the cytoskeleton, whereas GAPDH, EF1a, UBQ, and other proteins are involved in fundamental metabolism, which suggests that these genes are stably expressed in all cells and physiological states [44,45]. Later studies showed that expression of housekeeping genes is not necessarily consistent in different species under the same stress condition [46]. For example, in this study, TUB was stably expressed after high-temperature treatment of mushrooms, and Fb15 and UBQ5 were found to be the best reference genes for gene expression studies in rice endosperm under high temperature [47]. Regarding the selection of optimal reference genes of plants under high-temperature stress, the research conducted by Yan Haidong and others [48] revealed that eIF4A, TBP-1, and E2 could be used as reference genes for perennial ryegrass under high-temperature stress. 26S rRNA [49] and EF1a [50] have the most stable expression in wheat and maize, respectively, after high-temperature treatment. In addition, SmEF1a and SmUCP represent the most stably expressed reference genes of the eggplant heat-sensitive strain 05-1 and heat-resistant strain 05-4, respectively, at various time intervals during high-temperature stress [47]. To our knowledge, the selection and validation of reference genes under hightemperature stress in L. edodes have not been evaluated. In this study, expression levels as well as the stability of 10 candidate reference genes were measured in L. edodes using 3 different algorithms. The three algorithm yielded quiet different results with respect to the ranking of the 10 candidate reference genes, indicating the importance of using more than one software to achieve the best results. Due to the  different algorithms of the three analysis software, there was a slight difference in the ranking of candidate reference genes, but with the comprehensive analysis of TUB, it is still the best choice as a reference gene under high-temperature stress. In order to verify the result, two genes hsp100 and hyd1, which are related to fungal resistance to high-temperature stress, were selected. The two best candidate reference genes and  Figure 5: Validation of the selected reference genes. Relative expression of hsp100 and hyd1 in mycelia of L. edodes strains 18 (a; b) and 18N44 (c; d). In the strain 18, TUB and ACT were the most stable candidate reference genes, and UBC was the least stable candidate reference gene. In the strain 18N44, TUB and TUA were the most stable candidate reference genes, and GTP was the least stable candidate reference gene. Data are shown as mean ± standard deviation (n = 3). the worst candidate reference gene were used, respectively; the quantitative results of hsp100 and hyd1 were compared and verified. It was found that, when the two most stable candidate reference genes (TUB and ACT in strain 18; TUB and TUA in strain 18N44) were used, they have a similar expression trend, but while the least stable candidate reference gene (UBC in strain 18 and GTP in strain 18N44) was used, its expression of hsp100 or hyd1 had significantly higher than that of the two most stable internal reference gene, which would further prove the accuracy of TUB as a good reference gene for L. edodes under high-temperature stress. These results further implied that the use of an unstable reference gene could lead to the misinterpretation of gene expression results [51]. Therefore, choosing the best reference gene largely depends on the relevant factors such as species and tissue type, as well as experimental conditions. Different tissues and organs, growth and developmental stages, and biotic or abiotic stress as well as hormonal status will often, if not always, exert a certain influence on the specific expression of reference genes.

Conclusion
This study investigated the expression stability of reference genes of different L. edodes strains that have undergone hightemperature stress over different time periods. The 10 genes (TUB, TUA, GADPH, EF1, 18S, GTP, ACT, UBI, UBC, and