Identification and evaluation of reference genes for qRT-PCR studies in Lentinula edodes

Lentinula edodes (shiitake mushroom) is a common edible mushroom with a number of potential therapeutic and nutritional applications. It contains various medically important molecules, such as polysaccharides, terpenoids, sterols, and lipids, were contained in this mushroom. Quantitative real-time polymerase chain reaction (qRT-PCR) is a powerful tool to analyze the mechanisms underlying the biosynthetic pathways of these substances. qRT-PCR is used for accurate analyses of transcript levels owing to its rapidity, sensitivity, and reliability. However, its accuracy and reliability for the quantification of transcripts rely on the expression stability of the reference genes used for data normalization. To ensure the reliability of gene expression analyses using qRT-PCR in L. edodes molecular biology research, it is necessary to systematically evaluate reference genes. In the current study, ten potential reference genes were selected from L. edodes genomic data and their expression levels were measured by qRT-PCR using various samples. The expression stability of each candidate gene was analyzed by three commonly used software packages: geNorm, NormFinder, and BestKeeper. Base on the results, Rpl4 was the most stable reference gene across all experimental conditions, and Atu was the most stable gene among strains. 18S was found to be the best reference gene for different development stages, and Rpl4 was the most stably expressed gene under various nutrient conditions. The present work will contribute to qRT-PCR studies in L. edodes.


Introduction
The basidiomycete Lentinula edodes (Berkeley) Pegler (Lentinus edodes), an edible mushroom (shiitake), is the most economically important cultivated mushroom in East Asia and the second most popular in the world [1]. L. edodes is valued for its nutritional and medicinal properties, and has wide culinary and industrial applications [2]. Various compounds with therapeutic properties have been isolated from the mycelium and fruiting body in the last few decades, among which the polysaccharide lentinan has been studied intensively. Lentinan

Sampling and culture conditions
Xin 808, a common artificially cultivated strain in China and other Asia countries, was used. The fungus was grown on potato dextrose agar plates for 7 days at 25˚C. Three mycelial discs from the actively growing peripheral region were inoculated into 50 mL of the given medium. To collect samples reared under different nutrient conditions, two different media were used, i.e., carbon and nitrogen media. The carbon media consisted of the following components (g L -1 ): peptone, 2.5; KH 2 PO 4 , 0.1; MgSO 4 , 0.1; vitamin B1, 0.01; and three different carbon sources (glucose, sucrose or starch), 2.5. The nitrogen media consisted of the following components (g L -1 ): glucose, 2.5; KH 2 PO 4 , 0.1; MgSO 4 , 0.1; vitamin B1, 0.01; and two different nitrogen sources (yeast extract or glycine), 2.5. The cultures were statically incubated at 25˚C in the dark and collected after 30 days.
Four strains were included: three cultivated strains (Xin808, Xiang240 and Susheng69) and one wild variety (LSM-1). The vegetative mycelia of these strains were cultured in liquid potato dextrose medium. The culture conditions were the same as those described above.

Total RNA extraction and cDNA synthesis
The frozen samples were ground to a fine powder in liquid nitrogen, and 100 mg of the material was used for RNA isolation. Total RNA was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. The purity and quantity of samples were measured using a Nano spectrophotometer (ND-1000 Thermo Scientific, Waltham, MA, USA), and integrity was checked by agarose gel electrophoresis. First-strand cDNA was synthesized from 2.0 μg of total RNA using a reverse transcription kit (TaKaRa, Shiga, Japan). The synthesized cDNA was diluted 10 times with nuclease-free water and stored at -20˚C.
Since sequence information for L. edodes is limited, publicly available gene sequences from Gymnopus luxurians (http://genome.jgi.doe.gov/Gymlu1/Gymlu1.home.html), the most closely related species to L. edodes, were used to obtain gene sequences for L. edodes [25]. The specific reference gene sequences of L. edodes were obtained from the genome database of L. edodes (http://legdb.chenlianfu.com/), using G. luxurians homologus DNA sequences (CDS) as a queries.

Quantitative real-time PCR
Real-time amplification reactions were performed in 96 well plates using SYBR Green detection chemistry and run in triplicate on 96-wells plates with the iCycler iQ5 thermocycler (Bio-Rad, USA). Reactions were prepared in a total volume of 20 μL containing: 1μL of 10 fold diluted template, 0.5μL of each amplification primer (1μM), 5μL SsoFast™ Eva-Green1 Supermix (Bio-Rad, USA), and nuclease-free water to a final volume of 20 μL. Non-template controls were also included for each primer pair. Reactions were carried under the following conditions: 95˚C for 2 min, followed by 40 cycles of 95˚C for 15 s,63˚C for 20 s and 72˚C for 26 s. A melting curve analysis was generated by heating the amplicon from 60 to 90˚C to confirm that a single product was generated by each reaction. Thermal cycling, fluorescence data collection, and data analysis were performed using the iCycler iQ5 Thermo-cycler (Bio-Rad, Hercules, CA, USA) detection system according to the manufacturer's instructions.
Standard RT-PCR was performed for all primer pairs and amplification products of the expected size for each gene were confirmed by electrophoresis on a 2% agarose gel. The primer amplification efficiency was determined from a standard curve generated by serial dilutions of cDNA (10-fold each) for each gene in triplicate. Correlation coefficients (R 2 values) and amplification efficiencies (E) for each primer pair were calculated from the slope of the regression line by plotting mean Cq values against the log cDNA dilution factor in Microsoft Excel using the equation E (%) = (10 (1/-slope −1) Ã 100.

Data analysis
To facilitate the analysis of the expression stability of candidate reference genes, samples that differ with respect to developmental stage (Development), strain (Strain), and nutrient source (Nutrient), and all test samples (All). Three publicly available software tools, i.e., geNorm (version 3.5) [28], NormFinder (version 0.953) [29], and BestKeeper [30], were used to analyze expression variation among the ten candidate reference genes. Data were exported to Excel and converted to appropriate input files according to the software requirements.

Validation of reference genes analysis
To confirm the effectiveness of the selected reference genes for data normalization, expression levels of UDP-glucose pyrophosphorylase (UGPase), phosphoglucose isomerase (PGI) and phosphoglucomutase (PGM), which play important roles in the biosynthesis of polysaccharides, were examined in different developmental stages. According to the comprehensive ranking for samples obtained at various developmental stages, Btu and 18S were the most two stable genes and were selected as reference genes for the expression analysis. Additionally, the least stable reference gene, F-actin, was used to analyze the influence of the selection of an inappropriate reference gene. The qRT-PCR amplification conditions were the same as those described above. The relative expression levels of the three genes were standardized according to the formula Y (%) = 10 -(ΔCt/3) × 100, where ΔCt is the difference between the cycle threshold value of the target gene and the reference gene [31,32].

Selection and identification of candidate reference genes
A total of 10 candidate reference genes, including three genes (18S, 28S, and Pma) that have already been used as reference genes and seven new genes: Atu, Btu, Apt, F-actin, Rpl2, Rpl4, and Tsb were chosen in this study.
To determine the specificity of the primers used in this study, agarose gel electrophoresis and a melting curve analysis were performed following the qRT-PCR experiment. All primer pairs amplified a single PCR product of the expected size (Fig 1), without any nonspecific amplicons or primer dimers. The specificity of the amplicon was also confirmed by melting curves generated by qRT-PCR (Fig 1) and a sequencing analysis. A standard curve was generated using 10-fold serial dilutions of cDNA to calculate the gene-specific PCR efficiency. The amplification efficiencies for the 10 genes varied from 89% (F-actin) to 110% (Rpl4 and 18S), and the correlation coefficients ranged from 0.938 for Rpl4 to 0.996 for 28S. The T m ranged from 81.2˚C for Pma to 84.5˚C for Rpl2 (Table 1).

Expression profiles of reference genes
The expression stability of each candidate reference gene was evaluated by measuring Ct values for each individual sample pool. The expression profiles of these reference genes across all samples are shown in Fig 2. The Ct values for the ten reference genes ranged from 11.86 (18S in a maltose fermentation sample) to 34.89 (Apt in the primordium of L. edodes), showing a wide range of variation. The average Ct value for each reference gene was~24.11, with remarkable variation (15.00 for 18S to 30.73 for Apt). Greater than 84% of Ct values were between 16 and 30. The lowest gene expression variation among all the studied samples was observed for 28S (15.14-19.55), whereas the expression level of Apt was the most variable (21.88-33.25). The remaining eight reference genes exhibited intermediate variation (Fig 2). The wide range of Ct values for the ten reference genes indicated that none of the reference genes had a constant expression pattern in different samples. Therefore, in order to select a proper and stable reference gene to normalize gene expression data in L. edodes, it is necessary to evaluate the expression stabilities of these candidate reference genes under particular conditions.

Reference gene stability analysis
Three common algorithms, i.e., geNorm, NormFinder, and BestKeeper, were used to evaluate the ten reference genes. GeNorm (v3.5) yields M values by the stepwise exclusion of unstable genes, followed by the recalculation of M [28]. Stably expressed genes have M values of less than 1.5, and the most stable gene has the lowest M value [28]. The expression stability was further analyzed using Microsoft Excel-based NormFinder, by comparing the variation within and between user-defined sample groups [29]. NormFinder provides a stability value (SV) for each candidate gene, which is a direct measure of the expression variation. A lower SV value corresponds to higher gene expression [29]. BestKeeper finds the most stably expressed genes by estimating the coefficient of correlation (r) and BestKeeper Index (geometric mean) [30].
The most stable reference gene shows the lowest value r value and standard deviation (CV ±SD). Candidate genes that show a SD value of less than 1 are considered to be reliable references [19,[33][34][35]. The samples used in this study were divided into three subgroups to evaluate the expression stability of candidate genes using these three softwares. The first subgroup consisted of four divisions based on strains (Strain), the second subgroup consisted of three divisions based on the different developmental stages of Xin808 (Development), and the third subgroup consisted of five divisions based on the nutrient source of Xin808 (Nutrient). In addition, all samples were analyzed (All).

Reference genes expression among strains
The ranking of the ten reference genes from different strains obtained with the three algorithms indicated that Pma was the least stably expressed gene (S1 Fig). For most of the reference genes, the rankings obtained using NormFinder and geNorm were quite similar. However, greater variations was observed using BestKeeper, consistent with previous reports [15,17,36]. Each reference gene was evaluated with respect to stability using the three tests (Table 2), and Atu was the most stable gene, followed by Rpl4.

Reference gene expression among various developmental stages
The stability of the ten candidate reference genes among developmental stages was evaluated for all three samples. F-actin exhibited wide variation in Ct values (Fig 2) and was ranked the least stable gene according to all three algorithms (S2 Fig). 18S, was the most stable gene according to geNorm and was ranked second and fifth by NormFinder and BestKeeper, respectively. However, Pma, which was the least stable reference gene among strains, was identified as the fourth most stable gene among developmental stages. In all cases, 18S was the gene identified as the most suitable reference gene for the normalization of qRT-PCR data for

Reference gene expression for various nutrient conditions
The ranking of the reference genes for different nutrient conditions varied according to the three algorithms (S3 Fig). 28S was the least stable gene according to both geNorm and Norm-Finder, but was the most stable gene according to BestKeeper. Rpl4 was the best gene according to geNorm and NormFinder and was ranked fifth by BestKeeper. The three least stable genes were the same as those identified on the basis of expression profiles, i.e., F-actin, Pma, and Apt (Fig 2). Rpl4 was identified as the most stable gene, followed by Atu (Table 2). It may be noticed that the ranking of the ten candidate reference genes was quite different among sample groups.

Reference gene expression in all samples
To identify reference genes suitable for all the samples, the stability of the ten candidate reference genes was evaluated for all 11 samples and the results are presented in Fig 3. Pma and Factin were still identified as the two least stable genes by all three algorithms. Rpl4 and Btu, which were ranked as the top two genes by both geNorm and NorFinder, was identified as fifth and sixth most stable genes by BestKeeper. Using the stability ranking of each reference gene in all three tests, Rpl4 was identified as the most stable gene, followed by Btu and Rpl2 (Table 2). 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 [28,37]. To determine the optimal number of reference genes in each experimental condition, pairwise variation (V) was calculated using geNorm by applying a cutoff value of 0.150 [28]. A pairwise variation analysis showed that the optimal number of reference genes may be different for distinct samples. In the pairwise variation analysis for different developmental stage samples, the V2/3 value was less than 0.150, showing that only two reference genes is suffieient for normalizing gene expression in these samples (Fig 4). The pairwise variation among various nutrient sample conditions was V2/3 0.168, which is above the cutoff of 0.150, and the addition of next reference gene decrease this value to 0.130, implying that proper normalization required at least three reference genes. As a consequence, for nutrient samples, Rpl4, Atu and Tsb can be used as reference genes to normalize gene expression data. Similarly, the V5/6 value for various strains samples was 0.127, indicating that five reference genes are required to normalized gene expression strains. Accordingly, Atu, Rpl4, Btu, 18S, and Rpl2 can be used for data normalization in various strains. For all samples combined, the pairwise variation at V3/4 was 0.156, and with an additional gene it decreases to 0.128, indicating that at least four reference genes are needed for normalization. Therefore, Rpl4, Btu, Rpl2, and Tsb and can be used to normalize gene expression data for these samples.

Validation of reference genes
To validate the selection of candidate reference genes, the relative expression levels of three key genes (UGPase, PGI and PGM) involved in the biosynthesis of polysaccharides during various growth stages, including in the mycelium, primordium and fruiting body, were monitored. The results revealed that the three genes were differentially or specifically expressed (Fig 5).
The relative expression levels of UGPase and PGIwere higher in the fruiting-body than at the other two stages. No transcript of PGM was detected in the fruiting-body stage. Relative expression levels of PGM obtained by normalizing against "18S+Btu" or 18S alone indicated similar trends. However, normalization against F-actin led to the opposite conclusion, indicating that expression was higher in the fruiting-body stage.

Discussion
Among various technologies, qRT-PCR is considered the most appropriated method for gene expression analyses. However, several quality control measures can improve the reliability of results, including the validation of suitable reference genes for data normalization. In order to screen for suitable genes for data normalization in L. edodes, expression stabilities of ten house-keeping genes were analyzed using three software packages: geNorm, NormFinder, and Bestkeeper. The three statistical algorithms have been used to select and validate reference genes for qRT-PCR data normalization across a variety of tissues, organs, developmental stages, stress conditions in various plant species like tomato [20], pepper [21], bamboo [22], and fungi like (e.g.,Ganoderma lucidum) [18], and some specific conditions, e.g., in the rice blast fungus Magnaporthe oryzae [38] or plum under different postharvest treatments [24]. To our knowledge, the selection and validation of reference genes in L. edodes have not been evaluated.
In this study, expression levels as well as the stability of ten candidate reference genes were measured in L. edodes using three different algorithms. The candidate reference genes showed a relatively wide range of expression profiles under different conditions (Fig 2), confirming that no single gene had a constant expression profile. The three algorithm yielded quiet different results with respect to the ranking of the ten candidate reference genes, indicating the importance of using more than one software type to achieve the best results. Rpl4 was identified as the most stable gene for various nutrient conditions and all samples by NormFinder and geNorm, and was also recommended as the best reference gene for qRT-PCR data normalization in the entire sample in G. lucidum [18]. The NormFinder and geNorm methods ranked Atu as the first and second most stable genes among strains, respectively.18S was found to be the most stable gene according to geNorm, and NormFinder ranked 18S second, but Bestkeeper ranked it fifth. F-actin, a classic reference gene frequently used as an internal control for RT-PCR normalization in plants and animals, was identified as the least stable gene under nutrient, developmental variation and in all samples used in this study. Other than for the set of all samples, we can see that the top three reference genes were almost the same for all groups when determined by geNorm and NormFinder, and not by BestKeeper. This variation in the results may be explained by difference in the algorithms implemented by the three software packages.
It has been reported that the transcription levels of Apt and Pma remains constant throughout all stages of fruiting body development (mycelium, primordium, young fruiting body and mature fruiting body) in L. edodes [26], indicating that these two genes are stably expressed and can be used as an internal controls to study gene expression profiles for samples obtained at different developmental stage. However, Apt gene was ranked eighth for Development, Strain and Nutrient samples in this study (Table 2). This difference may be due to differences among strains used in the two studies. Pma, which has been used as an internal control to study fruiting body development in L. edodes [39], was ranked fifth by geNorm and NormFinder, and seventh by BestKeeper for different developmental samples ( Table 2), indicating that this gene is not suitable for developmental sample data normalization. These results emphasize the importance of validating reference genes for each experimental design.
In recent years, intensive research has focused on polysaccharide biosynthesis owing to its function as a host defense potentiator and thereby its important antivirus, antitumor, and immune function regulation roles among other functions. The polysaccharide content varied among developmental stages [40]. Hence, we analyzed the expression of three key genes involved in the polysaccharide biosynthetic pathway at different developmental stages. UGPase catalyzes the reversible production of uridine diphosphate glucose (UDPG) and pyrophosphate (PPi) from Glc-1-P and UTP, which is a key regulatory step for carbohydrate metabolism in Development sample (18S and Btu) and the least stable gene F-actin were used for normalization. For 18S and Btu geometric mean was calculated and used for normalization of expression. Error bars show standard deviation calculated from three biological replicates. The relative expression levels of the three genes were indicated as percentage to the expression of candidate reference genes. Different letters above histograms indicate statistical difference highlighted using ANOVA (p-value < 0.05). https://doi.org/10.1371/journal.pone.0190226.g005 Reference genes selections for Lentinula edodes many species because it is an universal precursor involved in biosynthesis of many carbohydrates, including sucrose, glucan, cellulose, hemicellulose, the carbohydrate moiety of glycolipids, and glycoproteins [41]. Higher expression levels of UGPase were detected in primordium stages than in other stages. PGM is involved in conversion of glucose-6-phosphate to glucose-1-phosphate, and there is a linear relationship between its activity and the production of EPS [42]. No transcript of PGM was observed in the fruiting-body stage, and its relative expression level in the primordium was higher than that in the fruiting-body stage Since PGI catalyzes the reversible isomerization of D-glucopyranose-6-phosphate and D-fructofuranose-6-phosphate, and the reaction that generates D-fructose-6-phosphate (F6P), it plays dual roles in EPS [43]. The transcript level of PGI was higher in the fruiting-body than in other developmental stages. The expression patterns of these three genes showed highly similar trends when the most stable reference genes 18S and Btu were used, either singly or in combination. However, a different conclusion, i.e., that the expression of PGI in the fruiting-body stage was higher than that in the other two stages, was obtained when the most unstable gene, F-actin, was used as a reference gene. These results further implied that the use of an unstable reference gene could lead to the misinterpretation of gene expression results.
In conclusion, the present work provides basic background information regarding reference gene selection for qRT-PCR studies in L. edodes. The study contributes to will make a contribution for further genomics and transcriptomics studies of this valuable medicinal fungus.