Selection and validation of reference genes for normalization of quantitative real-time reverse transcription PCR analysis in Poria cocos (Schw.) Wolf (Fuling)

Background Quantitative real-time reverse transcription PCR (qRT-PCR) requires a stable internal control to avoid misinterpretation of data or errors for gene expression normalization. However, there are still no validated reference genes for stable internal control in Poria cocos (Schw.) Wolf (Fuling). This study aims to validate the reference genes of P. cocos. Methods This study firstly collected the 14 candidate reference genes by BLASTP from the genome of P. cocos for qRT-PCR analysis to determine the expression levels of 14 housekeeping genes (GAPDH, MAPK, β-Act, RPB2, RPB1-1, RPB1-2, his3-1, his3-2, APT, SAMDC, RP, β-Tub, EIF, and CYP) under different temperatures and in response to different plant hormones (indole-3-acetic acid, abscisic acid, 6-benzylaminopurine, methyl jasmonate, and gibberellic acid), and the threshold cycle (Ct) values. The results were analyzed by four programs (i.e., geNorm, NormFinder, BestKeeper, and RefFinder) for evaluating the candidate reference genes. Results SAMDC, his3-2, RP, RPB2, and his3-1 were recommended as reference genes for treating P. cocos with indole-3-acetic acid, abscisic acid, 6-benzylaminopurine, methyl jasmonate, and gibberellic acid, respectively. Under different temperatures RPB2 was the most stable reference gene. CYP was the most stable gene for all 90 samples by RefFinder. Conclusion SAMDC, his3-2, RP, RPB2, and his3-1 were evaluated to be suitable reference genes for P. cocos following different treatments. RPB2 was the most stable reference gene under different temperatures and CYP was the most stable gene in the mycelia under all six evaluated conditions.


Background
Quantitative real-time reverse transcription PCR (qRT-PCR) is used for determining the abundance of mRNAs in molecular biology studies. Suitable reference genes are necessary to ensure accuracy and to avoid bias. Typically, reference genes are housekeeping genes necessary for cellular metabolism. The genes for cyclophilin (CYP), tubulin, ubiquitin, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), actin, 18S ribosomal RNA, 28S ribosomal RNA, and albumin are among the most frequently used reference genes [1].
However, the expression levels of reference genes may not be stable in different species [2], different tissues [3], or even identical cells under different culture conditions [4]. For example, the biosynthesis of triterpenes was induced by methyl jasmonate (MeJA) in Ganoderma lucidum (Leyss. ex Fr.) P. Karst (Lingzhi) [5,6]. However, the stability of fungal reference genes in the presence of plant hormones has not been properly evaluated by the gene expression levels of enzymes involved in the triterpene biosynthesis pathway.

Open Access
Chinese Medicine Little research has been conducted on reference genes in fungi. In Hemileia vastatrix Berk. and Br. (Toubaoxiujun), the cytochrome b, 40S ribosomal protein and Hv00099 genes have been selected as reference genes in vitro; however, the 40S ribosomal protein, GAPDH, and Hv00099 genes were the most stable genes in planta [7]. In Hypocrea jecorina Berk. and Br. (Hongherouzuojun), the gene encoding a GTPase was recommended as a reference gene [8]. Reference genes for qRT-PCR under different culture conditions and at different developmental stages in G. lucidum were reported [9].
Poria cocos (Schw.) Wolf (Fuling) is medicinal fungi and nutrition food widely distributes in East Asia, particularly in China, North America, Africa, and Australia [10,11]. Pharmaceutically active constituents extracted from P. cocos, including polysaccharides, triterpene derivatives, lanostane derivatives, and poricoic acid, exhibited anti-oxidant [12,13], anti-inflammatory [14], anti-tumor [15][16][17], anti-emetic [18], anti-nephritic [19], anti-rejection [20], diuretic [21], and anti-hyperglycemic activities [22]. The nematicidal activity of P. cocos was investigated and the active compounds were isolated [23]. Studies on the molecular biology of P. cocos were limited, including the basic molecular studies such as gene expression analysis and gene function identification [24]. qRT-PCR method was effective to detect the candidate genes involved in secondary metabolite biosynthesis. For example, the genes are most likely involved in the biosynthesis of pachymic acid in P. cocos was identified by qRT-PCR [25]; however, contigs and singletons were used instead of reference genes. The stability of potential internal control genes in P. cocos has not been evaluated.
This study aims to discover and obtain the stable reference genes of P. cocos for normalization of qRT-PCR analysis.

Sampling and culture conditions
The P. cocos strain CGMCC5.78 was purchased from the Institute of Microbiology, Chinese Academy of Sciences and was stored in the Institute of Medicinal Plant Development at −80 °C. We identified the strain using the DNA barcoding method with ITS2 primers. Ninety mycelial samples under different culture conditions were used in this study. Vegetative mycelia were cultured in two different media: potato dextrose agar medium (AOBOX, Beijing, China) and sucrose medium. The components of the sucrose medium were as follows: vitamin B1, 0.05 g/L; MgSO 4 •7H 2 O, 0.5 g/L; KH 2 PO 4 •H 2 O, 1 g/L; yeast extract, 2.5 g/L; peptone, 5 g/L; and sucrose, 35 g/L. The strain was maintained in potato dextrose agar medium. In the preculture stage, 40-mL sucrose medium was inoculated with mycelia and shaken (Thermo Fisher Scientific 491, Waltham, MA, USA) at 50 rpm in the dark in an incubator at 28 °C for 1 week. Subsequently, all of the mycelia were spread and were shaken at 120 rpm for an additional week in the dark at 28 °C. Finally, all cultures, including the culture broth, were incubated under various conditions (Table 1), including different concentration of hormones and different temperatures for 24 h.
The samples were arbitrarily allocated into six groups for analysis ( Table 1). The samples in groups A, B, C, D, and E were cultured in the media supplemented with different concentrations of indole-3-acetic acid (IAA; Sangon, Shanghai, China), abscisic acid (ABA; Sangon), 6-benzylaminopurine (6-BA;Sigma, St Louis, MO, USA), methyl jasmonate (MeJA; Sigma), and gibberellic acid (GA; Sangon), respectively. Group F comprised samples collected from cultures incubated at five different temperatures. The mycelia were collected by double gauze filters (CWBio, Beijing, China). Each experiment was performed in triplicate. A total of 90 samples were collected, and all of the samples were frozen in liquid nitrogen and stored at −80 °C.

Total RNA extraction, DNase treatment, and cDNA synthesis
The liquid nitrogen frozen samples were ground into fine powder by a mortar and pestle. The total RNA of each sample was extracted by the Polysaccharide and Polyphenol Total RNA Isolation Kit (spin column; BioTeke, Beijing, China) according to the manufacturer's instructions. The total RNA integrity and quality were confirmed by 1 % agarose gel electrophoresis by ethidium bromide staining. The RNA concentration was determined by a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). One microgram of total RNA of each sample was reverse transcribed by the Fast-Quant RT Kit (with gDNase; TIANGEN, Beijing, China) according to the manufacturer's protocol. All templates were diluted 30-fold for PCR and qRT-PCR.

Candidate gene selection, primer design, and validation
Based on previous studies [1,3,4] of reference genes determined in other species, 14 genes were evaluated in the present study, including multiple-copy genes. These genes include glyceraldehyde-3-phosphate dehydrogenase (GAPDH), mitogen-activated protein kinase (MAPK), beta actin (β-Act), RNA polymerase subu-nit2 (RPB2), RNA polymerase subunit1 (RPB1), histone 3 (his3), adenine phosphoribosyl transferase (APT), S-adenosyl methionine decarboxylase (SAMDC), ribosomal protein (RP), beta tubulin (β-Tub), eukaryotic translation initiation factor (EIF), and cyclophilin (CYP). The primer sequences, amplicon size and number of gene copies in the genome are summarized in Table 2. The candidate genes were selected from the P. cocos genome sequence database (SRA: PRJNA42921) by the BLASTP program (National Library of Medicine, USA) and a threshold E-value <1 × 10 −50 . Primer Premier 6.0 (PREMIER Biosoft, USA) and DNAMAN (LynnonBiosoft, USA) were used for primer design with the following criteria: an amplicon size ranging from 130 to 180 bp, an optimal T m of 53-55 °C, and a primer length from 18 to 22 bp. The primers were synthesized by Sangon Biotech (Shanghai, China). The specificity of each primer pair was measured by 2 % agarose gel electrophoresis following PCR (95 °C for 5 min; 35 cycles of 95 °C for 15 s and 60 °C for 1 min; 72 °C for 10 min) by the 90 cDNA sample mixture. Additionally, qRT-PCR was performed and the melting curve was determined for primers specific validation.

Real-time PCR performance and C t data collection
The expression level of each gene was determined in 96-well plates by an Applied Biosystems 7500 Real-Time PCR system (Life Technologies, Grand Island, NY, USA). Each reaction mixture contained 200 nM of each primer, 2 µL of the prepared cDNA template, 4.9-µL ddH 2 O, and 7.5-µL Ultra SYBR Mixture with ROX (CWBio, Beijing, China) in a final volume of 15 µL. The amplifications were performed by an initial denaturation step of 95 °C for 5 min, followed by 45 cycles of 95 °C for 15 s and 60 °C for 1 min. A temperature ramp step was added after 45 amplification cycles for specificity analysis (melting curve), with 95 °C for 15 s, 60 °C for 1 min, 95 °C for 15 s, and a final temperature of 60 °C for 15 s. There were three biological duplicate samples, and each biological duplicate sample was evaluated in triplicate.

Data analysis
The C t values from each reaction were used for analysis of the expression levels of all detected reference genes. The geNorm [26], NormFinder [27], BestKeeper [28], the Delta CT method [29] and the Web-based tool RefFinder [30] were used to determine the stability of the candidate reference genes. The default parameters of these software were applied.

Expression profile of candidate reference genes
The mean C t value was computed by three biological duplicates and three technical replicates for each independent experiment (the template generated from each condition of P. cocos was used in different independent experiment), and the three technical replicates were performed independently. A higher C t value indicates decreased transcription of the target gene. The average C t value of each candidate gene under conditions ranged from 22.45 ± 0.97 to 32.86 ± 0.86 cycles ( Table 3). The average C t value of six of the 14 genes was higher than 30.00. RPB1-2 and CYP demonstrated the lowest and highest relative expression levels, with average C t values of 31.21-33.21 and 22.37-23.91, respectively. The variation in the C t value was determined by the maximum and minimum C t values. The variation in the C t value of each candidate reference gene in all 90 samples ranged between 3.22 and 7.89. RPB1-1 exhibited the lowest variation in C t value followed by CYP (3.24). In contrast, EIF exhibited the highest variation in C t value (7.89).

Stability ranking of candidate reference genes
geNorm ranks the potential reference genes on the basis of their average pairwise variation in expression of one gene compared with each other gene of the set [26]. geNorm recommends 1.5 as the M-value cutoff. An M-value of less than 1.5 indicates stable expression, with the lowest M-value corresponding to the highest stability, and vice versa. Two reference genes were recommended for an ideal relative quantitative analysis. The M-values of candidate genes under different conditions generated by geNorm are listed in NormFinder is an Excel-based program for evaluating the expression stability of candidate reference genes based on the expression values, which enables estimation not only of the overall variation of the candidate normalization genes but also of variation between sample subgroups of the sample set [27]. NormFinder shows less sensitivity toward coregulation of the candidate normalization genes. A lower stability value indicates a higher stability. In group A, SAMDC was the most stable gene, with a stability value of 0.135, whereas his3-2 was the most unstable gene, with a stability value of 0.769. In group B, his3-2 exhibited the lowest stability value of 0.088, and β-Act exhibited the highest stability value of 1.586. Under different temperatures, his3-2 was the most stable, with a stability value of 0.069, and SAMDC was the least stable, with a stability value of 1.428. In group C, RP exhibited the best performance with a stability value of 0.183, and the expression level of RPB1-2 varied the most under different concentrations of 6-BA, with a stability value of 1.228. Following MeJA treatment, RPB2 exhibited the lowest variation, with a stability value of 0.106, and GAPDH exhibited the lowest stability value of 0.887. In group E, his3-1 was recommended as the reference gene for GA treatment, with a stability value of 0.244, and MAPK was the most unstable gene among the 14 genes, with a stability value of 1.119. When all of the samples were analyzed, CYP exhibited the lowest stability value of 1.320, whereas EIF exhibited the highest stability value at 13.240.
Gene expression stability was evaluated by BestKeeper using the standard deviation (SD), percentage covariance (CV), and correlation coefficient (r) [28]. BestKeeper can determine the best suited standards, out of 10 candidates, and combine them into an index. The candidate reference genes with SD >1 are considered unstable, and a higher  RefFinder analysis integrates four different methods (i.e., Delta CT, geNorm, NormFinder, and BestKeeper). The C t values were input into RefFinder directly, and the ranking of the four methods was calculated. Based on the rankings from each method, RefFinder assigns an appropriate weight to an individual gene and calculates the geometric mean of their weights for the overall final ranking [30]. The rankings of the candidate reference genes used in Delta CT were according to the repeatability of the gene expression differences among the samples. The results analyzed by RefFinder are summarized in Tables 5, 6 , 7, 8, 9, 10. In group A, SAMDC was recommended as the most stable reference gene. In group B, his3-2 exhibited the best performance, whereas his3-1 exhibited the best performance in group E. Under different treatment temperatures and different concentrations of MeJA, RPB2 maintained a stable expression level.    (Table 7). Following comprehensive analysis of all of the samples under the various conditions by Delta CT, geNorm, NormFinder, and BestKeeper, CYP was recommended as the reference gene. The results obtained using these different methods were not identical. In group C, RP was recommended as the most stable gene by all of these methods, whereas in group F, RPB2 was recommended as the reference gene by Delta CT, geNorm and RefFinder. However, in the remaining groups, Delta CT, NormFinder and RefFinder recommended the same gene as the reference gene; SAMDC, his3-2, RPB2 and his3-1 in groups A, B, D, and E, respectively. Following comprehensive analysis of all of the samples under the various conditions, CYP was recommended as the reference gene by Delta CT, Best-Keeper, NormFinder, and RefFinder, although not with geNorm. According to the above-mentioned results, Ref-Finder was likely the most comprehensive and scientific of these methods.

Evaluation of the combination of reference genes
Pairwise variation (V) determines the optimal number of control genes for normalization and proposes 0.15 as a cutoff value [26]. If the V n /V n+1 value is less than 0.15, the suitable gene number for normalization is n. Additional control genes were not necessary in the six groups except for group C (i.e., the 6-BA treatment group),as indicated by V 2 /V 3 values below 0.15 [26]. Three reference genes were recommended for group C, as indicated by a V 3 /V 4 value of 0.115, which is consistent with the M-value ranking for this group.

Discussion
Validation of the stability of candidate reference genes under different experimental conditions [31], with different tissues [32,33], at different stages, and in different species [34] is necessary. In the present study, EIF was the most unstable gene in P. cocos; however, EIF1 and EIF3 were recommended as reference genes in Ammopiptanthus mongolicus (Maxim. ex Kom.) S.H. Cheng [35]. In contrast, CYP was the most stable gene in leaves of Deschampsia antarctica É. Desv. [36] under three abiotic stresses (salt, cold, and PEG treatment), whereas the EF-1α gene was recommended for roots. In banana fruit, the expression levels of two widely used reference genes, actin and GAPDH, were not stable [34].
The candidate reference gene rankings for the individual groups evaluated in this study may differ slightly from the ranking for all samples because, under specific circumstances, more accurate rankings would be established.
Moreover, most of the M-values of the 14 genes were less than 1.5 except for SAMDC, EIF, and APT, indicating that most of the candidate reference genes were stable. As one of the least stable genes, the instability of APT has been reported in papaya under six experimental conditions [37]. It was contradictory that CYP was the best overall reference gene but did not exhibit the best performance in any single group. CYP was the most stable reference gene using Delta CT, NormFinder, and Best-Keeper but not geNorm (Table 11). In addition, CYP was the third-most stable reference gene by geNorm. Moreover, CYP frequently ranked among the top five reference genes (Tables 5, 6,7,8,9,10), particularly under 6-BA treatment, in which CYP exhibited the highest average M-value when using geNorm for analysis. In group C, CYP ranked firmly as the second-most stable reference gene. In contrast, the ranking of other candidate genes in the six groups varied greatly. A similar phenomenon has been observed in Ammopiptanthus mongolicus [35]. EIF1 and EIF3 were selected as reference genes across all of the samples, whereas these two genes were the most stable only under drought stress among the four evaluated abiotic stresses. Following acibenzolar-S-methyl treatment, the combination of CYP and eIF4B was most suitable as an internal control in Eucalyptus L'Hér. In addition to P. cocos and Eucalyptus [38], CYP has been selected as an internal control for several animal cells. In human peripheral blood, CYP was a more suitable housekeeping gene than β-Act and GAPDH [39]. CYP was also recommended as one of the reference genes for neurons of the central nervous system [40] and in atopic human bronchial epithelial cells [41]. Moreover, CYP was considered to be an RNA normalization control in rats [42].
NormFinder, BestKeeper and geNorm are widely used for selection of reference genes, although the results generated by the different methods may be slightly different [43,44]. Our results displayed the same tendency as those of previous studies [26][27][28][29][30]. Moreover, the validity of the results might be related to the materials used or even to potential experimental errors. The importance of systematic evaluation before candidate genes are used as reference genes, especially under different conditions were observed in the study.

Conclusion
SAMDC, his3-2, RP, RPB2, and his3-1 were evaluated to be suitable reference genes for P. Cocos following different treatments. RPB2 was the most stable reference gene under different temperatures and CYP was the most stable gene in the mycelia under all six evaluated conditions.