Selection and Evaluation of Potential Reference Genes for Quantitative Real-Time PCR in Agaricus blazei Based on Transcriptome Sequencing Data

Quantitative real-time PCR (qRT-PCR) is widely used to detect gene expression due to its high sensitivity, high throughput, and convenience. The accurate choice of reference genes is required for normalization of gene expression in qRT-PCR analysis. In order to identify the optimal candidates for gene expression analysis using qRT-PCR in Agaricus blazei, we studied the potential reference genes in this economically important edible fungus. In this study, transcriptome datasets were used as source for identification of candidate reference genes. And 27 potential reference genes including 21 newly stable genes, three classical housekeeping genes, and homologous genes of three ideal reference genes in Volvariella volvacea, were screened based on transcriptome datasets of A. blazei and previous studies. The expression stability of these genes was investigated by qRT-PCR analysis and further evaluated by four software packages, geNorm, NormFinder, BestKeeper, and RefFinder. Among these candidates, α-TUB (Tubulin alpha) and Cox5a (COX5A subunit VA of cytochrome c oxidase) were revealed as the most stable in fruit body, and suitable for 5 different developmental stages. α-TUB and ATP3 (ATP3 gamma subunit of the F1 sector of mitochondrial F1F0 ATP synthase) showed the most stable expression in stipe tissues and, Uqcrc (core subunit of the ubiquinolcytochrome c reductase complex) and PUP3 (20S proteasome subunit beta 3) performed well in pileus tissues during the process of A. blazei development, while GAPDH (glyceraldehyde-3-phosphate dehydrogenase) was among the least stable genes in all sample sets. Finally, the Ableln3 (homology of eln3 gene of Coprinus cinereus) was adopted to validate the reliability of these stable and unstable reference genes, indicating that the use of unsuitable reference genes as internal controls could change the target gene’s expression pattern. This study can provide guidance for choosing reference genes for analyzing the expression pattern of target genes and facilitate the functional genomic investigation on fruit body formation and development, as well as stipe elongation and pileus expansion in A. blazei.


Introduction
Quantitative real-time PCR (qRT-PCR) has been commonly used as a powerful technique for rapid measuring gene expression level and validate transcriptomic data [1,2]. However, accurate and reliable qRT-PCR analysis requires suitable internal control genes to normalize gene expression [3,4]. Ideally, the expression of reference genes is unaffected by various biological or experimental conditions, such as different tissues, different developmental stages, or experimental treat-ment [1,5]. A handful of traditional reference genes, for instance GAPDH, β-actin, and 18S rRNA, has been prevalently used in the literature [6][7][8]. Recent studies, however, showed that the stability of traditional internal control genes varied considerably under different experimental conditions in various species [9][10][11][12]. The application of an improper reference gene for normalization in qRT-PCR analysis could bring about misinterpretation of expression data, with consequent incorrect results [13,14]. Therefore, the systematic stability analysis of reference genes should be evaluated throughout all development stages and tissues. For accurate and reliable qRT-PCR analysis, it is necessary to select appropriate internal control genes to normalize gene expression.
With the rapid development of high-throughput sequencing technology, RNA-sequencing has been widely employed to analyze transcriptome in various species [15][16][17]. Benefiting from this technology, a large number of differentially expressed genes (DEG) and presumptive function genes associated with reproductive and development were identified [16,18,19]. Additionally, transcriptome sequence data were also applied as a source of robust normalization genes [1,20,21]. In the study of Narsai et al. [1], 151 genes with consistent expression level under all conditions were identified and, 12 of these genes were then selected and validated with qRT-PCR for different tissues and under stress. In Oxytropis ochrocephala, 12 genes were chosen as candidate reference genes on the basis of the abiotic transcriptome datasets by Zhuang et al. [21].
Agaricus blazei Murrill, popularly known as Cogumelo do Sol in Brazil and in China as Ji Song Rong, is one of the most valuable edible and culinary-medicinal mushroom species [22][23][24]. Primarily, a housekeeping gene GAPDH has been employed as a reference gene for normalization in A. blazei [24]. However, the stability of this gene has not been estimated in this fungus. Recently, genome and transcriptome of A. blazei were investigated to identify genes involved with fruit body development and biosynthesis pathway of polysaccharide and benzaldehyde [25,26]. The transcriptome datasets could be used as a potential source for selection of appropriate candidate reference genes in A. blazei.
In this study, 21 candidate reference genes were selected based on the transcriptome datasets of A. blazei. In addition, three classical housekeeping genes, GAPDH, α-tubulin, and β-tubulin, and homologous genes of three ideal reference genes (SPRYp, Ras, and Vsp26) in Volvariella volvacea were also included. The expression stabilities of twenty-seven candidates were carefully assessed using qRT-PCR analysis and four popular programs. Additionally, Coprinus cinereus eln3 homolog (Ableln3) was selected as target gene to evaluate the effectiveness of selected reference genes.

Materials and Methods
2.1. Strain, Incubation Condition, and Sample Harvest. A. blazei heterkayon strain JA was obtained from Fujian Academy of Agricultural Sciences and then cultured on potato dextrose agar (PDA) at 24°C. Spawn of strain JA was cultured on straw-manure compost [26]. For RNA extraction, samples of mycelia (MY) and fruit bodies at different developmental stages (primordia (PR), harvested stage (HA), and pileus opening stage (OP)) were collected according to the method described by Lu et al. [25,26]. Fruit bodies at button stage (BU) and tissue samples (stipe and pileus from fruit bodies at button stage, harvested stage, and pileus opening stage) were also harvested, respectively, and frozen in liquid nitrogen immediately. Three independent biological replicates were conducted for each sample.

Screening of Stably Expressed A. blazei Genes and
Traditional Internal Control Genes. Transcriptome data across four different developmental stages (Table 1) [25,26] were calculated by an FPKM formula. To assess expression stability, we calculated the following indices of FPKM values [21], including standard deviation (SD), mean expression value (MV), and coefficient of variation value (CV, dividing SD by MV) for each gene in accordance with the method by de Jonge et al. [27]. Genes with lower CV values were thought to be more stable, and the cut-off CV value was set as less than 0.3 [1,21,28]. Additionally, DEG profiles of that differentially expressed (Padj < 0:05 and |log 2FoldChange | ≥1) between any two sets of the 8 examined samples, and genes predicted as hypothetical genes or genes with unknown function were filtered [29]. The candidate reference genes with high similarity (identify ≥ 90) with other fungi were chosen for further study.
2.3. RNA Extraction and cDNA Synthesis. Total RNA was isolated with RNAprep Pure Plant Kit (Polysaccharides and Polyphenolics-rich) (Tiangen, Shanghai, China) in accordance with manufacturer's instructions. Additionally, NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) were employed to evaluate RNA purity and concentration. Synthesis of cDNA was conducted using TransScript All-in-One First-Strand cDNA Sythesis SuperMix for qPCR kit (TransGen Biotech, Beijing, China) according to the manufacturer's protocol.

Primer Design and qRT-PCR.
With the aim of estimation of gene expression level, primer premier 5.0 [30] was applied to design primer pairs (Table 1) for qRT-PCR analysis. And qRT-PCR reaction was performed using TB GreenTM

Samples
Accession number  Samples  Accession number  Samples  Accession number   MY1  GSM3814130  HA3  GSM3814147  OP2  GSM3814155  MY2  GSM3814131  HA_pileus1  GSM3814148  OP3  GSM3814156  MY3 GSM3814132 Premix Ex TaqTM (Tli RNaseH Plus) (TaKaRa, Dalian, China) with CFX Connect real-time system (Bio-Rad, USA). And then, 2% agarose gel electrophoresis was used to confirm the specific product of each primer. Standard curves of five points and five-fold dilution series from pooled cDNA were generated [31] and adopted to calculate PCR efficiency (E). To verify the reliability of the selected reference genes, we evaluated the relative expression pattern of Ableln3 gene (a homolog of C. cinereus eln3 involved in fruiting body morphogenesis [32] (Ableln3F: GGAGAATGTTTCCGCAGGT, Ableln3R: CCACACTACGAGGACACGCT). And the 2 −△△Cq method was employed to calculate the relative expression [33].

Data Processing and Gene Expression Stability Analysis.
The experiment dataset was divided into four subsets. Subset A included mycelia and samples from four different developmental stages of fruit body. Subset B was similar to subset A but absence of mycelia. Subset C and subset D were composed of stipe and pileus from fruit bodies at button stage, harvested stage, and pileus opening stage, respectively. Additionally, we analyzed the overall stability of 27 genes among total samples. And four software programs, geNorm, NormFinder, Bestkeeper, and RefFinder, were applied to assess the stability of the candidate genes. For geNorm [34] and NormFinder analysis [35], the arithmetic mean Cq (quantification cycle) value obtained from three biological and three technological replicates [29,36] were transformed to relative quantities using the equation Q = 2 −ΔCq (ΔCq = each corresponding Cq value − the lowest Cq value) [37]. We then imported relative quantities into geNorm and calculated the expression stability (M value) of each gene. M value less than 1.5 indicated that the gene can be selected for a reference gene, and the gene with the lowest M value was considered as the most stable gene [36,37]. And then, pairwise variation (Vn/Vn + 1) value was employed to analyze the optimal number of reference genes, and the values <0.15 (the recommended cut-off threshold) indicated that the number of the most suitable reference genes was n [36,37]. Otherwise, an additional reference gene was required for normalization. Similar to geNorm, the stability value (SV) of genes was calculated by NormFinder, and the lower SV means the more stable of the gene. BestKeeper analysis, in which geometric mean of Cq values were used directly, was performed to investigate the gene expression stability by calculating CV, SD, and correlation coefficient between Cq value and Bestkeeper index (calculated as the geometric mean of Cq values of candidates) [38]. Finally, RefFinder (https://omictools.com/ reffinder-tool) was employed to make a combination of the outcome from three methods mentioned above and to generate an integrated ranking of the most suitable reference genes in different sample sets [37,39].

Identification of Candidate Reference Genes.
In present study, a total of 27 genes (Table 2) were obtained. Our gene panel contained (i) a group of 21 genes with predicted function (identify ≥ 90%) on the basis of their stable expression patterns across mycelia and the different developmental stages of fruit body as well as across tissues samples (stipe and pileus from fruit body at harvested stage and pileus opening stage); (ii) housekeeping genes (such as GAPDH, α-tubulin, and β-tubulin2); and (iii) homologous genes of three ideal reference genes (SPRYp, Ras, and Vsp26) in V. volvacea [29].

Quantitative Real-Time PCR Analysis of Candidate
Reference Genes in A. blazei. The difference in expression levels between 27 candidate reference genes was evaluated based on qRT-PCR. The specificity of qRT-PCR amplification for each candidate was verified by melting curve and 2% agarose gel electrophoresis analysis. The results showed that all of the primer sets generated single band with expected size, ranging from 89 to 331 bp ( Figure S1). The PCR products were then further confirmed by sequencing and the specificity were also validated by the present of single peak in melting curve analysis ( Figure S1). The amplification efficiency (E) of each primer pair ranged from 90.0% for Vps26 to 104.3% for ATPDH ( Table 2). The correlation coefficient (R 2 ) varied from 0.986 (ATPDH) to 1.000 (Uqcr8). Thus, the primer sets were acceptable for further qRT-PCR analysis.
In general, the exploratory analysis of Cq values ( Figure 1) displayed the expression levels of the reference genes and offered us a glimpse into gene stability. The average Cq values of 27 candidates, calculated across all samples, revealed a maximum of 25:72 ± 0:91 (RPT6) and minimum of 19:94 ± 0:92 (ATP2) (Figure 2) for lowest and highest expression levels, which showed the expression level of RPT6 was approximately 50-fold lower than ATP2. The relative expression variation (the maximum Cq value minus the minimum Cq value) for each candidate showed the lowest variation for NdufA2 indicating that it was the most stably expressed in total samples and could be chosen as the optional reference gene. Meanwhile STT3 displayed greatest variation and should be avoided as a reference gene. Therefore, 26 genes except STT3 were selected for further analysis.

Expression Stability of Candidate Genes in A. blazei.
In view of the complex surroundings of edible mushroom, the stability of reference genes under different conditions needs to be explored systematically. And four different software programs NormFinder, geNorm, BestKeeper, and RefFinder were adopted to estimate the expression stability of candidate reference genes [5].
3.3.1. GeNorm Analysis. M values, negatively correlated with gene stability, were calculated using geNorm software. GeNorm analysis showed that all reference genes possessed an M value less than 1.5 indicating that all genes performed well under individual sample sets. For different developmental stages, Cox5a and α-TUB were the most stable genes among the 26 candidates. Cox5a and ATP1 were the two best genes within the developmental stages of the fruit body and total samples. For the pileus and stipe of fruit bodies at different developmental stages, the most stable genes were Uqcr8 and PUP3, ATP1, and ATP3, respectively, whereas the traditional reference genes GAPDH and β-TUB2 in all sample sets with 3 BioMed Research International  (Figure 2).

NormFinder
Analysis. The NormFinder program was also applied to determine the stability of reference genes, and the results were shown in Table 3. The outcomes from geNorm and NormFinder analyses were somehow similar. The twelve most invariable reference genes were almost similar except set C (stipe of fruit body), though the ranking orders were different in the results obtained by the two approaches ( Figure 2 and Table 3). The most appropriate gene obtained by NormFinder analysis in all of the sample sets was almost in the top nine appropriate genes generated by geNorm analysis. Additionally, the most variable genes were the same in both geNorm and NormFinder analyses in all sample sets. genes in set C and set D (stipe and pileus of fruit body at different developmental stages, respectively) ( Table 4). And 24 genes had SD value smaller than 1.0 when total samples were taken into consideration. Furthermore, genes with a higher R and a significant P value (P value < 0.05) were identified as the most stable candidates [42,43]. Based on this principle, ATP3, Copb1, FKS1, and ATP3 were the most stable candidate reference genes in 5 stages, different developmental stages of fruit body, stipe, and pileus of fruit body.

3.3.4.
RefFinder. According to RefFinder, a web-based analysis tool, the most stable genes were generated in the five developmental stages (Figure 3(a)), different developmental stages of fruit body (Figure 3(b)), stipe (Figure 3(c)), pileus (Figure 3(d)), and total samples (Figure 3(e)). α-TUB, Cox5a, ATP7, and ATP2 were found to be four of the top five stable genes in the different developmental stage of fruit body and five developmental stages. α-TUB was also the top three appropriate gene in remaining sample sets except set 3 (pileus) in which this gene was ranked at 18, while our results showed that GAPDH was the least stable gene in all sample sets.

Determination of the Optimal Number of Reference
Genes. The optimal number of the candidate reference genes was also determined using geNorm based on pairwise variation (Vn/Vn + 1) in each sample set. As shown in Figure 4, all   BioMed Research International the Vn/Vn + 1 value in all sample sets was lower than 0.15. Thus, 2 reference genes would be sufficient for qRT-PCR data normalization in A. blazei.

Reference Genes Validation.
In order to evaluate the reliability of candidate reference genes, the expression levels of Ableln3 were investigated across 5 developmental stages and different developmental stages of fruit body using qRT-PCR analysis with α-TUB and Cox5a. Meanwhile, the expression patterns of Ableln3 were also detected in different tissues of fruit body such as stipe (α-TUB and ATP3) and pileus (Uqcrc and PUP3). Moreover, we evaluated the expression of Ableln3 following normalization with one of the least stable reference gene (GAPDH) for a comparative analysis. It was observed that the expression pattern of Ableln3 was similar by using α-TUB and Cox5α across different   BioMed Research International developmental stages of fruit body and 5 developmental stages ( Figure 5). However, the relative expression of target gene Ableln3 was underestimated in fruit body at OP stage, when normalized using GAPDH. Additionally, similar expression patterns were obtained using α-TUB and ATP3 in stipe tissues and Uqcrc and PUP3 in pileus tissues, respectively. Whereas, GAPDH used for normalization, the target gene's expression patterns were significantly different from the above patterns in pileus tissues.

Discussion
Due to its high sensitivity, specificity, and high-throughput feature, qRT-PCR has become a commonly and powerful tool for detecting gene expression [1,2,44,45]; however, it requires suitable reference genes for data normalization which has been considered as one of the vital factors affecting the reliability and accuracy of the quantitative gene expression analysis [2,43]. Validation of proper candidate reference genes is mandatory to yield reliable results. Therefore, numerous studies have focused on reference gene stability to identify the most suitable reference genes [37,[46][47][48].
Most of them, however, investigate the stability of traditional housekeeping genes such as GAPDH, actin, 18S rRNA, tubu-lins, and elongation factor 1 alpha. Genes encoding GAPDH had been used as internal control genes for A. blazei in previous studies [24], whereas there is an absence of analysis of reference genes based on developmental stage and tissue type. Hence, a systematic exploration and evaluation of candidate reference genes is required for study on development in A. blazei.
To achieve an accurate conclusion, more than two programs were recommended for stability evaluation of reference genes because different statistical methods based on distinct calculating principles might generate different conclusion from the same data [49]. For example, only two programs (geNorm and NormFinder, the most widely used algorithms) were adapted to study the stability of gene expression by Kanakachari et al. in Solanum melongena; however, the results yield from the two programs showed some differences [4,50]. Thus, the use of at least three different algorithms was required to achieve reliable results [49,51]. In the present study, twenty-seven candidate reference genes were chosen based on A. blazei transcriptome data and previous studies and systematically analyzed using geNorm, NormFinder, BestKeeper, and RefFinder during different developmental stages and in different tissue samples of fruit body, to determine suitable candidate reference genes in A. blazei. The results obtained from geNorm and Norm-Finder were somehow similar, but there were some differences from the results generated by Bestkeeper. In BestKeeper, Cdc4 ranked fourth in set A, while it occupied a medium position by geNorm and NormFinder. And Cox5a was the most stable candidate gene identified by geNorm and NormFinder in set B, whereas it appeared at the eleventh place in the BestKeeper results. The results were similar to many previous studies, and the small discrepancies were probably due to the different principles in the three statistical algorithms [41,51,52]. RefFinder, an integrative statistical program, has been widely employed to calculate the overall stability of candidate reference genes and determine suitable reference genes for diverse species [48,52,53]. Based on the comprehensive RefFinder analysis, α-TUB and Cox5a, ATP7 and α-TUB, α-TUB and ATP3, and Uqcrc and PUP3 were indicated as the most appropriate internal control genes for normalization in fruit body, 5 developmental stages, stipe tissues, and pileus tissues, respectively. The integrated results yielded from different statistical programs could lead to better accuracy for each gene, suggesting that at least three algorithms should be used for stability evaluation of candidate reference genes.
In this study, the expression stability of commonly used reference genes including GAPDH, α-TUB, and β-TUB and three most stable expressed genes (Ras, Vsp26, and SPRYp) in V. volvacea were tested. Amongst these selected traditional genes, α-TUB was identified as top three stable reference genes in 5 developmental stages, different developmental stage of fruit bodies, stipe tissues, and total samples, but it was ranked as one of the least stable gene in pileus. The expression stability of β-TUB, GAPDH, SPRYp, and Ras were ranked late, and GAPDH was found to be unsuitable candidate reference genes in pileus tissues. The results reported here highlighted the fact that the traditional reference gene may not always display invariable expression; therefore, it is necessary to evaluate their stability prior to qRT-PCR analysis in different species.
Recently, increasing studies suggested that more than one reference gene should be used as internal control genes because the use of more than one gene could avoid the erroneous information which might be provided by using only single one [48,54,55]. Based on the geNorm pairwise variation (Vn/Vn + 1), the minimum number of candidates required for normalization was two across all sample sets, which was similar with previous studies by Li et al. [11] and       Figure 4: Determination of the optimal number of reference genes for normalization using geNorm. Pairwise variation (Vn/n + 1) less than 0.15 indicated an additional reference gene was not required.
Yang et al. [48]. Based on validation of target gene's expression, α-TUB and Cox5a were recommended as appropriate reference genes for qRT-PCR normalization in fruit body, and they were also suitable for 5 developmental stages. ATP3/α-TUB and Uqcrc/PUP3 could be used as stable reference genes for normalization in stipe tissues and pileus tissues, respectively.
In summary, this study illustrated that transcriptome sequence data were useful as a source of reference gene screening. Based on our transcriptome databases and previous studies, we selected 27 reference gene candidates and evaluated the expression stability of these genes in different developmental stages and tissue samples (stipe and pileus) of fruit body in A. blazei. α-TUB and Cox5a were identified as the optimal candidates for fruit body and appeared at the top three positions in 5 different developmental stages. For stipe tissues, ATP3 and α-TUB were the most suitable reference genes. And the use of Uqcrc and PUP3 was sufficient to provide accurate relative expression quantification in pileus tissues. These results can be helpful to improve the accuracy of gene expression using qRT-PCR analysis and facilitate the functional genomic investigation on fruit body formation and development, as well as stipe elongation and pileus expansion in A. blazei.

Data Availability
The transcript raw Illumina sequencing data of all samples were submitted to NCBI Gene Expression Omnibus (GEO) and the accession numbers have been described in Table 1.  Figure S1: validation of primer specificity by 2% Agarose gel electrophoresis and melting curves. (a) Agarose gel electrophoresis displayed the amplification of a single PCR product of the expected size for 27 candidate reference genes and Ableln3, M indicated DL2000 bp DNA marker. (b) The melting curves for 27 candidate reference genes and Ableln3 showed single peaks. (Supplementary Materials)