Selection and Validation of Reference Genes for qRT-PCR in Lentinula edodes under Different Experimental Conditions

Lentinula edodes is the most consumed mushroom in Asia due to its nutritional and medicinal values, and the optimal reference gene is crucial for normalization of its gene expression analysis. Here, the expression stability of 18 candidate reference genes (CRGs) in L. edodes was analyzed by three statistical algorithms (geNorm, NormFinder and BestKeeper) under different stresses (heat, cadmium excess and Trichoderma atroviride infection), different substrates (straw, sawdust and corn stalk) and different development stages (mycelia, primordia and fruit bodies). Among the 18 CRGs, 28S, Actin and α-tub exhibited the highest expression stability in L. edodes under all conditions, while GPD, SPRYP and MSF showed the least stable expression. The best reference gene in different conditions was different. The pairwise variation values showed that two genes would be sufficient for accurate normalization under different conditions of L. edodes. This study will contribute to more accurate estimation of the gene relative expression levels under different conditions using the optimal reference gene in qRT-PCR (quantitative reverse transcription polymerase chain reaction) analysis.


Introduction
With the advantages of high sensitivity and repeatability as well as dynamic quantification range, quantitative reverse transcription polymerase chain reaction (qRT-PCR) is widely used to analyze and validate the expression of target genes [1][2][3]. Nevertheless, the results of qRT-PCR were usually not accurate because of errors caused by many factors such as the efficiency of complementary DNA (cDNA) synthesis and amplification, gene expression normalization, and so on [4,5]. In order to minimize those errors, one or more genes which are usually stably expressed in any experimental conditions are applied in qRT-PCR analysis [5,6]. However, the increasing evidence demonstrated that the expression levels of many tested internal genes in plants and animals as well as fungi were relatively constant only in specific cells or experimental conditions [7][8][9]. The application of unstable internal reference genes will result in significant variations or differences in the qRT-PCR results [10]. Hence, it is necessary to identify and analyze the stability of the reference genes under different conditions to guarantee the accuracy and reliability of qRT-PCR results.
For hundreds of years, Lentinula edodes has been used as decoctions and essences or as alternative medicine in many Asian countries such as China, Korea, Japan, etc. [11]. In addition, L. edode is the most widely cultivated edible fungi in the world because of its nutritional properties and effective bioconversion of agricultural wastes [12,13]. As a white-rot fungus, the genes in L. edodes were identified to

Strain, Culture Conditions and Sample Collection
In treatment with different carbon sources, L. edodes strain W1 (ACCC50926) was cultured separately at 25 • C on sawdust medium [21], straw medium (33.3 g straw, 20 g wheat bran, 2 g gypsum and 20 g agar in 1 L sterilized water) and corn stalk medium (33.3 g corn stalk, 20 g wheat bran, 2 g gypsum and 20 g agar in 1 L sterilized water). For abiotic and biotic stress treatment group, mycelia of L. edodes strain W1 cultured on sawdust medium were treated under 40 • C heat stress for 24 h, and mycelia growing on the PDA (potato dextrose agar) medium were subjected to Trichoderma atroviride infection for 24h. Meanwhile, mycelia cultured on PDA medium 5 mg/mL Cd(NO 3 ) 2 were also collected. Meanwhile, primordia and fruiting bodies of L. edodes strain W1 were also collected according to Xiang's report [27]. After collection, all samples were immediately frozen in liquid nitrogen and stored at −80 • C for RNA extraction.

RNA Isolation and cDNA Synthesis
Total RNA was extracted from all samples using the RNAiso Plus Kit (TaKaRa, Dalian, China) according to the manufacturer's instruction. The RNA integrity and quantity were checked by 1% agarose gel electrophoresis and the NanoDrop DS-11 Spectrophotometer (Applied Denovix, Madison, WI, USA). The RNA samples with the OD 260 /OD 280 value between 1.8 and 2.2 were reverse transcribed into cDNA using HiScript ® II Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech, Nanjing, China) according to the instructions of the manufacturer.

Selection of Candidate Reference Genes and Primer Design
Based on previous studies [26,27], sequences of eighteen genes, including seven house-keeping genes and eleven novel candidate genes, were downloaded from the genome database of L. edodes [16] and used for primer design. The detailed information of these genes was displayed in Table 1. Primers were designed using Primer Premier 5.0, and the length of PCR amplification products varied from 100 to 253 base pairs (Table 1).

qRT-PCR, Amplification Efficiency and Data Analysis
The qRT-PCR reaction was performed using CFX Connect Real-Time system (Bio-Rad, Hercules, California, USA). The PCR amplification mixture contained 1 µL of cDNA, 5 µL of AceQ ® qPCR SYBR ® Green Master Mix (Vazyme Biotech, Nanjing, China), 2.5 µL of 10 nM primer mixture and 1.5 µL ddH 2 O. The PCR reaction was performed with an initial denaturation for 5 min at 95 • C; 40 cycles of 30 s at 95 • C, annealing at 60 • C for 30 s; an extension at 72 • C for 20 s. The melting curve (60-95 • C) was used to determine the specificity of every qRT-PCR reaction. All qRT-PCR reactions were performed in biological and technical triplicates. The standard curves were generated by qRT-PCR detection using a 10-fold dilution gradient of the first cDNA (1, 10 -1 , 10 -2 , 10 -3 , 10 -4 ) as templates with each repeated three times. Amplification efficiencies (E) and correlation coefficients (R 2 values) were checked from standard curves. The qRT-PCR Ct values of eighteen genes were converted to relative copy numbers using the standard curves, and used to carry out geNorm and NormFinder and Bestkeeper analysis.

Validation of Reference Genes
To validate the reliability and accuracy of the selected reference genes for data normalization, the expression levels of seven genes were analyzed in different experimental conditions (Table S1). The expression data of the target genes were normalized with the two most stable genes and the least stable reference genes. The qRT-PCR amplification conditions were the same as those described above. The relative expression levels of all target genes were calculated by the 2 −∆∆Ct method.

Analysis of Candidate Reference Genes under Different Conditions Based on RNA-Seq Data
To analyze the expression stability of the eighteen CRGs, the RPKM (reads per kilobase per million mapped reads) values were gained from RNA-seq data in each experiment involved in different L. edodes strains, biotic/abiotic stresses and substrates. The RNA-seq data of 24 h Trichoderma atroviride stress in L. edodes strains YS3334 and YS55 and different carbon substrates for W1 were obtained from our lab studies (unpublished data). The expression data from heat (S606 and YS3357) and light (L135, ACCC 50903) treatments were obtained from previous studies [18,20]. The expression variation of the eighteen CRGs was determined using a threshold of an absolute log 2 (fold change between treatment and control RPKM values) ≤ 1.

Expression Profile of Candidate Reference Genes in L. edodes
The cycle threshold (Ct) values were used to normalize for the standard curves generated with ten 10-fold serial cDNA dilutions by qRT-PCR. The formula E = 10 -1/slope -1 was used to calculate the amplification efficiency (E) varied from 90.6% for 28S to 102.3% for MSF1 (Table 1). All correlation coefficients (R 2 ) of eighteen genes were greater than 0.99, indicating a reliable linear relationship of the respective Ct values with the log values of the initial gene copy numbers. The mean Ct values of the eighteen candidate genes in all samples under different conditions varied from 15.25 to 30.84 (Figure 2A), and ten of the candidate genes had a Ct mean values lower than 25, including Actin, GPD, α-tub, β-tub, EF, 28S and UBI, demonstrating that the expression levels of they were very high in all samples. EF showed considerably higher expression than the other candidate reference genes (low Ct values), while SPRYP exhibited the lowest expression (highest Ct values) in all samples. Relative variation analysis of the Ct values (the maximum Ct value minus the minimum Ct value) in all candidate reference genes documented the greatest variation in SPRYP while the smallest variation in GPD of L. edodes.

Expression Profile of Candidate Reference Genes in L. edodes
The cycle threshold (Ct) values were used to normalize for the standard curves generated with ten 10-fold serial cDNA dilutions by qRT-PCR. The formula E = 10 -1/slope -1 was used to calculate the amplification efficiency (E) varied from 90.6% for 28S to 102.3% for MSF1 (Table 1). All correlation coefficients (R 2 ) of eighteen genes were greater than 0.99, indicating a reliable linear relationship of the respective Ct values with the log values of the initial gene copy numbers. The mean Ct values of the eighteen candidate genes in all samples under different conditions varied from 15.25 to 30.84 (Figure 2A), and ten of the candidate genes had a Ct mean values lower than 25, including Actin, GPD, α-tub, β-tub, EF, 28S and UBI, demonstrating that the expression levels of they were very high in all samples. EF showed considerably higher expression than the other candidate reference genes (low Ct values), while SPRYP exhibited the lowest expression (highest Ct values) in all samples. Relative variation analysis of the Ct values (the maximum Ct value minus the minimum Ct value) in all candidate reference genes documented the greatest variation in SPRYP while the smallest variation in GPD of L. edodes.

Evaluation of Expression Stability of the Eighteen Candidate Reference Genes in L. edodes
According to qRT-PCR Ct values, three statistical algorithms geNorm and NormFinder and BestKeeper were used to evaluate the stability of the eighteen selected CRGs. All selected genes were used for geNorm and NormFinder analyses, while, for BestKeeper analysis, only the top ten stable genes obtained with the two previous algorithms were used. GeNorm analysis of candidate reference genes. The expression stability measure (M) of a CRG was calculated as previously reported [36]. The smaller the M values are, the higher the expression stability will be. The values of the eighteen CRGs varied from 0.758 to 0.086 ( Figure 2 and Table S2). From the perspective of different development stages, β-tub and UBI were the most stable genes with the smallest value of 0.249, followed by Actin and RPA12, in contrast to the least stable genes of GPD and CYPL with M values of 0.723 and 0.646, respectively ( Figure 2B). Under the conditions of different substrates, α-tub and β-tub were most stable with the smallest value of 0.086, while CYPL and MSF were less stably expressed ( Figure 2C). In terms of different stresses, α-tub and Actin were the most stably expressed with an M value of 0.159, whereas SPRYP and MSF were less stably expressed

Evaluation of Expression Stability of the Eighteen Candidate Reference Genes in L. edodes
According to qRT-PCR Ct values, three statistical algorithms geNorm and NormFinder and BestKeeper were used to evaluate the stability of the eighteen selected CRGs. All selected genes were used for geNorm and NormFinder analyses, while, for BestKeeper analysis, only the top ten stable genes obtained with the two previous algorithms were used. GeNorm analysis of candidate reference genes. The expression stability measure (M) of a CRG was calculated as previously reported [36]. The smaller the M values are, the higher the expression stability will be. The values of the eighteen CRGs varied from 0.758 to 0.086 ( Figure 2 and Table S2). From the perspective of different development stages, β-tub and UBI were the most stable genes with the smallest value of 0.249, followed by Actin and RPA12, in contrast to the least stable genes of GPD and CYPL with M values of 0.723 and 0.646, respectively ( Figure 2B). Under the conditions of different substrates, α-tub and β-tub were most stable with the smallest value of 0.086, while CYPL and MSF were less stably expressed ( Figure 2C). In terms of different stresses, α-tub and Actin were the most stably expressed with an M value of 0.159, whereas SPRYP and MSF were less stably expressed ( Figure 2D). All the integrated data suggested that 28S and 18S were the most stably expressed genes, while GPD and SPRYP were the least stably expressed genes ( Figure 2E). In addition, the pairwise variation values (Vn/Vn+1) calculated by geNorm were used to determine the optional number of reference genes for accurate normalization. As reported previously, an additional reference gene would be unnecessary with the pairwise variation value less than 0.15 [37]. As shown in Figure 2F, all pairwise variation values were smaller than 0.15, demonstrating that two genes are sufficient for accurate normalization under different conditions of L. edodes. The optimal combination of reference genes was UBI with β-tub for different development stages, β-tub with α-tub for different substrates and α-tub with Actin for different stresses, respectively. Nevertheless, 28S with 18S was the optimal combination of reference genes under all conditions (Table S2).
NormFinder analysis of candidate reference genes. The CRGs were ranked by NormFinder based on the intra and inter group expression variations [36]. A notable difference was observed between the results of NormFinder and geNorm analyses about the most and least stably expressed genes. In NormFinder analysis, UBI and RPA12 were predicated as the most stable genes for different substrates and different development stages, respectively, while α-tub was defined as the most stable gene under different stress and total conditions (Table 2). Meanwhile, ranked in the second and third positions were RPA12 and 18S for different substrates, α-tub and Actin for different development stage, 18S and PPCI for different stresses, and 28S with β-tub for total conditions. The least stable genes were CYPL for different substrates, SPRYP for different stresses and GPD for different development stages and total conditions. BestKeeper analysis of candidate reference genes. The top ten CRGs ranked by the sum of stability values in geNorm and NormFinder analyses were selected for BestKeeper analysis. The expression variability of the top ten CRGs was measured by the standard deviation (SD) and correlation coefficient (R). The lower the values of SD and CV were, the higher the stability of CRGs would be. The CRGs with a SD value larger than 1 were considered as unstably expressed genes [38]. As shown in Table 3, RPA12 and 28S were defined as the most stable genes for different substrates, while 18S, PPCI and UBC with a low R value were identified as unstably expressed genes. In terms of different development stages, the most stable gene was Actin, followed by RPA12 and UBC ranked in the second and third positions, respectively. Meanwhile, 28S was identified as the most stably expressed gene under different stresses, followed by PPCI and RPA12 in the second and third positions, respectively. Nevertheless, 18S was eliminated due to its high p-value (0.08866) and low R value (0.6380). The CRGs with the low p-values (0.001) and a high correlation coefficient (>0.93) were defined as the most stably expressed genes, including four genes (RPA12, PPCI, 28S and Actin). Based on the results of the three statistical analyses, Actin and RPA12 were defined as the two most stable genes in different stages, while 28S with RPA12 and α-tub with Actin were defined as the most stably expressed for different substrates and different stresses, respectively. In summary, 28S could be treated as the reference gene under all different conditions.

Validation of Selected Reference Genes in L. edodes
To validate the reliability of the optimal and least reference genes for qRT-PCR analysis under different conditions, the relative expression levels of seven genes were calculated.
From the perspectives of growth and development stages, it was documented that the relative expression levels of three tested genes exhibited a significantly difference using different reference genes ( Figure 3A). In comparison to the mycelia stage, LeDnaJ11 normalized by GPD was obviously upregulated under primordium and fruiting body stage, whereas there was no difference found in the analysis with Actin and 28S as the reference genes. Meanwhile, we found that Lelcc1 and Lelcc11 in primordium stage were significantly downregulated when normalized by Actin and 28S, but no obvious difference was detected with the reference gene GPD.

Evaluation of the Expression Stability of Candidate Reference Genes in Different L. edodes Strains
To estimate the stability of the selected CRGs in different L. edodes strains, the RPKM values were gained from RNA-seq databases [20,37]. Genes with the RPKM ratio between any two treatments displaying a less than 2-fold change between treatment and control (non-treated) samples were considered as stably expressed. Under heat stress, for L. edodes strain YS3357, VPS28 and PPCI displayed an obvious upregulation, in contrast to the downregulation of CYPL and PP2A; for L. edodes strain S606, four genes (UBI, PPCI, VPS28 and CYPL) were upregulated, in contrast to the downregulation of SPRYP (Figure 4 and Table S3). From the perspective of T. atroviride infection, sixteen of eighteen CRGs were stably expressed except for UBI (downregulated in YS55) and MSF (upregulated in YS3334). During the coloring process in the presence or absence of light, only one gene RPA12 showed upregulation in L. edodes strain L135. When the mycelia of L. edodes strain W-1 mycelia were cultured in the CYM liquid medium containing different carbon sources (glucose, cellulose or cellulose and sodium lignosulphonate), all of the CRGs were stably expressed ( Figure 4 and Table S3). In terms of different substrates, an obvious difference was found in qRT-PCR analysis normalized by different reference genes ( Figure 3B). The relative expression level of Lelcc9 in straw medium exhibited a similar upregulation trend using different reference genes, but the upregulation range with GPD was significantly higher than those of Actin and 28S. In addition, it was observed that the relative expression levels of LeDnaJ11 and LePOD-Mn normalized by GPD showed an obvious upregulation in straw medium, whereas no difference was found in the analysis with Actin and 28S as the reference genes.
Under different stress conditions, a prominent difference was detected in qRT-PCR analysis normalized by different reference genes ( Figure 3C). It was found that the relative expression levels of LeDna98 normalized by GPD showed an obvious upregulation when subjected to T. atroviride infection, yet no difference was observed in the analysis with Actin and 28S as the reference genes. In addition, the relative expression levels of LeDnaJ07 and LeHsp98 in heat stress showed a similar change trend using different reference genes, where the upregulation range with GPD was higher than those of Actin and 28S.

Evaluation of the Expression Stability of Candidate Reference Genes in Different L. edodes Strains
To estimate the stability of the selected CRGs in different L. edodes strains, the RPKM values were gained from RNA-seq databases [20,37]. Genes with the RPKM ratio between any two treatments displaying a less than 2-fold change between treatment and control (non-treated) samples were considered as stably expressed. Under heat stress, for L. edodes strain YS3357, VPS28 and PPCI displayed an obvious upregulation, in contrast to the downregulation of CYPL and PP2A; for L. edodes strain S606, four genes (UBI, PPCI, VPS28 and CYPL) were upregulated, in contrast to the downregulation of SPRYP (Figure 4 and Table S3). From the perspective of T. atroviride infection, sixteen of eighteen CRGs were stably expressed except for UBI (downregulated in YS55) and MSF (upregulated in YS3334). During the coloring process in the presence or absence of light, only one gene RPA12 showed upregulation in L. edodes strain L135. When the mycelia of L. edodes strain W-1 mycelia were cultured in the CYM liquid medium containing different carbon sources (glucose, cellulose or cellulose and sodium lignosulphonate), all of the CRGs were stably expressed ( Figure 4 and Table S3).

Discussion
Abiotic and biotic stresses and cultivation materials are the major limiting factors for the yield and quality traits of L. edodes. To solve these problems, we need to understand the molecular mechanism of L. edodes in response to stresses and the breeding of stress-resistant varieties. qRT-PCR can help us to understand this mechanism by quantifying the gene expression, and a stable internal gene reference gene is a prerequisite to ensure the reality of qRT-PCR results. Recently, an increasing number of studies have been performed on reference genes in organisms. For instance, several housekeeping genes were found to be unstably expressed under abiotic and biotic stresses [39,40], suggesting the necessity to identify reliable reference genes for qRT-PCR analysis. From the perspective of macro basidiomycetes, the stabilities of novel and traditional candidate genes have been estimated undertaken in G. lucidum, V. volvacea, P. ostreatus and L. edodes under different

Discussion
Abiotic and biotic stresses and cultivation materials are the major limiting factors for the yield and quality traits of L. edodes. To solve these problems, we need to understand the molecular mechanism of L. edodes in response to stresses and the breeding of stress-resistant varieties. qRT-PCR can help us to understand this mechanism by quantifying the gene expression, and a stable internal gene reference gene is a prerequisite to ensure the reality of qRT-PCR results. Recently, an increasing number of studies have been performed on reference genes in organisms. For instance, several house-keeping genes were found to be unstably expressed under abiotic and biotic stresses [39,40], suggesting the necessity to identify reliable reference genes for qRT-PCR analysis. From the perspective of macro basidiomycetes, the stabilities of novel and traditional candidate genes have been estimated undertaken in G. lucidum, V. volvacea, P. ostreatus and L. edodes under different development stages, strains, nutrient conditions and abiotic stresses [22][23][24][25][26][27]40]. These reports suggested that the use of reference genes should be based on specific stress conditions and sample tissues as well as strains.
In the present study, a comprehensive analysis of eighteen CRGs, including traditional and novel reference genes from previous studies, was performed under different stresses (40 • C heat, Cd 2+ excess and T. atroviride infection), different substrates (straw, sawdust and corn stalk) and different development stages (mycelia, primordia and fruit bodies). The melting curves as well as R 2 and E values exhibited high specificity and amplification efficiency for all eighteen CRGs (Figure 1 and Table 1).
Under heat and Cd 2+ as well as T. atroviride stress conditions, α-tub was demonstrated as the best reference gene by the results from the three statistical analyses, which was consistent with previous studies reporting that α-tub was stably expressed under many kinds of abiotic stresses in fungi (V. volvacea and L. edodes), animals and plants [23,26,[41][42][43]. PPCI and Actin were identified as the second most stable genes by the three statistical methods. Previous reports documented Actin as a stably expressed gene in the response of L. edodes to 37 • C stress and in tobacco and Jute subjected to abiotic stresses [26,44,45] and CYC3 (PPCI) as a reliable reference gene for Morchella sp. in temperature stress [24]. The pairwise variation V2/3 value was below 0.15, demonstrating that the two genes were suitable for normalization in gene expression analysis. Therefore, α-tub with Actin or PPCI was suggested as the pair of reference genes for L. edodes in response to different kinds of stresses.
In different substrates, 28S and RPA12 were determined by the three analysis methods as the most suitable primer pair of reference genes. The analysis of geNorm, NormFinder and BestKeeper ranked 28S in the fourth, third and second position, respectively. RPA12 was estimated as the most stable gene by BestKeeper analysis and ranked in the second position by geNorm analysis. For the entomopathogenic fungus Beauveria bassiana, 28S was confirmed as the most stably expressed gene under a number of nutritional and stress conditions [46]. However, RPA12 encoding DNA-directed RNA polymerase I subunit has rarely been reported as an optimal reference gene in previous studies. Therefore, 28S was recommended as the best internal control gene for the response of L. edodes to different substrates.
For different development stages, the selection of reference genes has been reported in different varieties of edible fungi. In V. volvacea, Ras and SPRYP were reported as the most stable genes in heterokaryon H1521 [40], while L-asp and MSF were considered as stably expressed genes in homokaryon PYd15 [23], which were different from the present study in that Ras and SPRYP as well as MSF were found to be unstably expressed in all samples. Additionally, CYC3 encoding peptidyl-prolyl cis-trans isomerase exhibited a stable expression in ten Morchella species [24]. In L. edodes strain Xin 808, 18S and Rpl4 (ribosomal protein L4) were confirmed as the stably expressed genes [27], while in, L. edodes strain, W-1 Actin and RPA12 were defined as the stably expressed genes in different development stages. Furthermore, Actin was suggested as the most relevant reference genes during berry development [47], and ACT1 was reported as the most suitable reference gene during rice seed development [48]. All these results indicate that Actin can be used as a reliable reference gene during L. edodes development.

Conclusions
Overall, the three different statistical algorithms identified α-tub as the optimal internal control gene for many kinds of stresses and different L. edodes strains, while Actin and 28S as the suitable reference genes for gene expression normalization in different development stages and substrates. This study contributes to accurate analysis of differential expression changes in L. edodes under different conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/9/647/s1, Figure S1: Amplification fragments of eighteen candidate reference genes by agarose gel electrophoresis, Table S1: Delineations of genes used for validation analysis of selected reference genes, TableS2: Stability analysis of eighteen CRGs was calculated by geNorm, TableS3