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

Abstract

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 ubiquinol-cytochrome 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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. 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 treatment [1, 5]. A handful of traditional reference genes, for instance GAPDH, β-actin, and 18S rRNA, has been prevalently used in the literature [6–8]. Recent studies, however, showed that the stability of traditional internal control genes varied considerably under different experimental conditions in various species [9–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–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–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.

2. 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.

2.2. 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 ($P\text{adj}<0.05$ and $\mid \log 2\text{FoldChange}\mid \ge 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 ($\text{identify}\ge 90$) with other fungi were chosen for further study.

Table 1

Transcriptome data of A. blazei.

Samples	Accession number	Samples	Accession number	Samples	Accession number
MY1	GSM3814130	HA3	GSM3814147	OP2	GSM3814155
MY2	GSM3814131	HA_pileus1	GSM3814148	OP3	GSM3814156
MY3	GSM3814132	HA_pileus2	GSM3814149	OP_pileus1	GSM3814157
PR1	GSM3814133	HA_pileus3	GSM3814150	OP_pileus2	GSM3814158
PR2	GSM3814134	HA_stipe1	GSM3814151	OP_pileus3	GSM3814159
PR3	GSM3814135	HA_stipe2	GSM3814152	OP_stipe1	GSM3814160
HA1	GSM3814145	HA_stipe3	GSM3814153	OP_stipe2	GSM3814161
HA2	GSM3814146	OP1	GSM3814154	OP_stipe3	GSM3814162

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.

2.4. 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 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^{-△△\text{Cq}}$ method was employed to calculate the relative expression [33].

2.5. 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^{-\Delta \text{Cq}}$ ($\Delta \text{Cq}=\text{each}\text{corresponding}\text{Cq}\text{value}-\text{the}\text{lowest}\text{Cq}\text{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 ($\text{Vn}/\text{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].

3. Results

3.1. 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 ($\text{identify}\ge 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].

Table 2

Selected candidate reference genes, primers, and amplification characteristics derived from qRT-PCR analysis.

Gene	A. blazei gene ID	Description	Primer sequence (5${}^{\prime }$-3${}^{\prime }$)	Product size (bp)	Amplification efficiency (%)	$R^2$
Uqcr8	A08809	Ubiquinol-cytochrome c reductase complex subunit 8	F:CACCCAATATGCCATTTCTCCR:AATCCCGTATCCGACACCAA	136	99.2	1.000
Vps26	A09130	Vacuolar protein sorting 26	F:GGGAGAAGATGCTCGGAAACR:CCGCCAACTCTTGCGATAA	230	90.0	0.998
Gnb2	A04205	Guanine nucleotide-binding protein beta subunit 2	F: ACAAAGTCCACGCCATCCCTR: GAGACCACCGCAAGCAACAT	89	92.9	0.998
Copb1	A08816	Coatomer subunit beta	F:TTATTGGGAAGTCTGTCCGAAGTR: TGAACGACAGGTCGGTATGAG	190	93.7	0.998
Uqcrc	A01805	Core subunit of the ubiquinol-cytochrome c reductase complex	F: ACAACCTCGCCAACTCATTCAR:TGCTCATTCTTGCCCATTCTTT	141	101.4	0.997
ATP7	A05732	ATP7 subunit D of the stator stalk of mitochondrial F1F0 ATP synthase	F: TCGCTCATTATCGCACTATCCTR:GGCACTGTCCATTTACCCTTCTT	319	93.7	0.999
NdufS7	A01315	NdufS7 ubiquinone oxidoreductase subunit	F: CCTTGAAGTCCGCAGTTTGTAATGR: TTTGTTTGTGGAGACGAGGTGAAT	100	97.6	0.998
PUP3	A04280	20S proteasome subunit beta 3	F:GAAGAGCGTGAGATTGAACCAGAR:TCTGTCGCAGCAATAAAGGGTT	152	95.2	0.996
PsmB6	A05770	Proteasome subunit beta type-6	F: CGTCCACGACCGCATTTACTR:AGGTCCACCCGACATTTGC	112	96.8	0.999
Ras	A03682	Ras protein	F:GGCAACAAGTGCGACAAATCAR:CCACAGCGAAAGCAACAGAAA	251	103.4	0.997
CPR1	A01938	Probable proteasome subunit beta type-4	F: CCTGCTGCTGCCTCGTCTTR:CTCTTGAGTGTTTCTATGCCTTCC	260	92.7	0.996
CS	A02079	Citrate synthase	F:GCTTGTCGGTTCAGCCCTTTCTR:TCTTCCTTAGAACGGCGTGGC	230	97.4	0.997
SPRYp	A04044	Endosomal SPRY domain protein	F: ACCCTATCCCAGTTCCTTTCCAR:CTTCACCTCCCAGTAGTAGACTTCATT	198	103.0	0.998
GAPDH	A07487	Glyceraldehyde-3-phosphate dehydrogenase	F: TTCGGCTCATTTGAAGGGTGR:TTGGTTGTGCAAGAAGCGTTAG	140	94.0	0.998
ATP2	A06017	ATP2 beta subunit of the F1 sector of mitochondrial F1F0 ATP synthase	F:CCGAATACTTCCGTGACGAGGR:TTGGTAACCGACAGCAGAGGG	128	99.0	0.999
ATP1	A07930	ATP1 alpha subunit of the F1 sector of mitochondrial F1F0 ATP synthase	F: TCCCTCCTTCTCCGTCGTCCR:TCTTAGTCTGGGCGGCTGAA	300	96.6	0.999
Cdc4	A08640	Myosin regulatory light chain cdc4	F: CAGGCAGAAGTCGCTGAGATTGTR:GAGGATTGTGCGGACAAAGGAT	300	96.3	0.992
FKS1	A09827	1, 3-beta-glucan synthase	F: ACCATTCACGCCGATTACATTR:GCATAGCATTACGCCAACGA	187	104.0	0.996
Cox5a	A03337	COX5A subunit VA of cytochrome c oxidase	F: TCTCGTCCGTTGTGCGTCTGR:TCGCGGTAGCGTAGTCGTTGA	206	93.2	0.996
ATPDH	A06239	ATP-dependent helicase	F: TCAAATCTGGCGTTGTCCCR:CCGATTCCTTGATGTTCATTACTC	211	104.3	0.986
α-TUB	A05905	Tubulin alpha	F:CATTTCCCACTTGCCACTTTCR:CAACAGCAGCATTGACATCCTTA	202	100.3	0.999
NdufA2	A00266	NdufA2, NADH ubiquinone oxidoreductase 10.5kD subunit	F: GCGAACTCCGTTTCTTCTTCTCCCR:AGTCTGGGTTGTGCTGTTTGATGG	102	96.2	0.994
RPT6	A01698	26S proteasome regulatory complex ATPase RPT6	F: GTTGTTTGTTATGGCTCGTGAAR:AAGCAACTCCAACATCGTCCT	136	96.3	0.993
ATP3	A08690	ATP3 gamma subunit of the F1 sector of mitochondrial F1F0 ATP synthase	F: ACAAGTTCGTTTCTGCCATCTCCR:TGGTAGCATCATCCTCGTTCTCAT	114	97.8	0.999
NDUFB9	A09944	NDUFB9 NADH-ubiquinone oxidoreductase complex I	F:TACCTCACAAACTCGCTTGACTGGR:CGGTTACATCACGGTTCTTCTCG	100	92.8	0.994
β-TUB2	A00521	Tubulin beta-2	F:CGGGCAAGATGTAGAGCACR:TCCAGGATAAAAGGCGAAG	231	100.8	0.996
STT3	A07729	Oligosaccharyl transferase STT3 subunit	F: TGTAGCCAACACTGCCTCCTR:TGCCAAAAGCACACAAATG	331	98.5	0.997

F and R indicated forward primer and reverse primer, respectively; $R^2$ referred to coefficient of determination.

3.2. 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\pm 0.91$ (RPT6) and minimum of $19.94\pm 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.

[figure omitted; refer to PDF]

[figures omitted; refer to PDF]

3.3. 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 high $M$ value were identified as unstable reference gene among 26 candidates (Figure 2).

3.3.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.

Table 3

Ranking of candidate reference genes yielded from NormFinder programs.

Rank	5 stage		Fruit body		Stipe		Pileus		Total
	Gene	Stability	Gene	Stability	Gene	Stability	Gene	Stability	Gene	Stability
1	α-TUB	0.044	α-TUB	0.049	ATP7	0.048	Uqcrc	0.075	ATP3	0.076
2	Cox5a	0.062	Cox5a	0.067	NdufS7	0.052	ATP7	0.083	α-TUB	0.078
3	Vps26	0.068	Vps26	0.074	Uqcr8	0.070	Copb1	0.088	Cox5a	0.082
4	ATP3	0.069	ATP3	0.078	Copb1	0.073	FKS1	0.088	ATP7	0.086
5	Uqcrc	0.080	ATP7	0.079	FKS1	0.076	NDUFB9	0.098	Vps26	0.087
6	ATP2	0.080	ATP2	0.088	Vps26	0.083	Uqcr8	0.099	ATP2	0.087
7	ATP7	0.081	NdufS7	0.089	α-TUB	0.084	ATP3	0.100	FKS1	0.089
8	NdufS7	0.083	Uqcrc	0.089	ATP3	0.093	ATP2	0.100	ATP1	0.089
9	ATP1	0.088	Copb1	0.094	RPT6	0.094	RPT6	0.109	RPT6	0.091
10	RPT6	0.091	RPT6	0.097	Cdc4	0.099	Cox5a	0.110	Uqcrc	0.092
11	PUP3	0.092	ATP1	0.099	Ras	0.101	PUP3	0.119	Copb1	0.097
12	CS	0.093	Uqcr8	0.102	Uqcrc	0.103	ATP1	0.120	NdufS7	0.097
13	Cdc4	0.094	CS	0.102	NDUFB9	0.110	NdufA2	0.123	Cdc4	0.097
14	FKS1	0.097	PUP3	0.103	PUP3	0.110	ATPDH	0.134	NDUFB9	0.103
15	Copb1	0.099	Cdc4	0.104	Cox5a	0.113	α-TUB	0.140	PUP3	0.105
16	Uqcr8	0.102	NDUFB9	0.106	ATP1	0.114	Cdc4	0.143	Uqcr8	0.111
17	NDUFB9	0.104	FKS1	0.109	CS	0.121	CPR1	0.145	CS	0.112
18	ATPDH	0.117	ATPDH	0.131	ATPDH	0.123	Vps26	0.149	ATPDH	0.117
19	Gnb2	0.131	Gnb2	0.148	CPR1	0.127	CS	0.157	CPR1	0.121
20	CPR1	0.139	CPR1	0.156	SPRYp	0.134	PsmB6	0.168	Gnb2	0.133
21	SPRYp	0.158	SPRYp	0.174	PsmB6	0.140	Ras	0.173	PsmB6	0.150
22	PsmB6	0.162	PsmB6	0.182	ATP2	0.142	NdufS7	0.173	Ras	0.151
23	β-TUB2	0.175	β-TUB2	0.196	Gnb2	0.148	Gnb2	0.212	SPRYp	0.167
24	Ras	0.188	Ras	0.211	GAPDH	0.213	β-TUB2	0.225	β-TUB2	0.187
25	NdufA2	0.223	NdufA2	0.212	β-TUB2	0.300	SPRYp	0.235	NdufA2	0.194
26	GAPDH	0.260	GAPDH	0.282	NdufA2	0.321	GAPDH	0.329	GAPDH	0.229

3.3.3. BestKeeper Analysis

BestKeeper program was then further employed to investigate the stabilities of candidate reference genes according to the standard deviation (SD) of the average Cq values, coefficient of variance (CV), coefficient correlation ($R$), and $P$ value. Genes with $\text{SD}\text{value}>1.0$ were considered to be unstable and unsuitable for gene expression normalization [5, 10, 40, 41]. BestKeeper analysis displayed that 19 genes in set A. (5 different developmental stages) had SD value smaller than 1.0, 25 genes in set B. (fruit body at different developmental stages), and 26 genes and 25 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.

Table 4

Expression stability of candidate reference genes obtained from BestKeeper analysis.

Rank	5 stages		Fruit body		Stipe		Pileus		Total
	Gene	SD [±CP]	Gene	SD [± CP]	Gene	SD [± CP]	Gene	SD [± CP]	Gene	SD [± CP]
1	ATP3	0.962	Copb1	0.627	FKS1	0.769	ATP3	0.368	ATP7	0.837
2	Vps26	0.788	Cdc4	0.564	NdufS7	0.647	Cdc4	0.418	ATP3	0.882
3	Uqcr8	0.978	ATP3	0.613	ATP7	0.656	α-TUB	0.409	Uqcrc	0.963
4	Cdc4	0.859	NdufS7	0.521	ATP3	0.631	Uqcrc	0.158	Copb1	0.906
5	α-TUB	0.981	ATP7	0.552	Uqcrc	0.668	Cox5a	0.401	Vps26	0.676
6	CS	0.962	α-TUB	0.534	Cdc4	0.693	SPRYp	0.527	Cdc4	0.843
7	NdufS7	0.811	NDUFB9	0.660	α-TUB	0.571	CS	0.370	Cox5a	0.828
8	ATP2	0.911	Vps26	0.387	Ras	0.731	CPR1	0.301	NdufS7	0.760
9	FKS1	0.909	Uqcr8	0.567	ATPDH	0.673	β-TUB2	0.527	FKS1	0.819
10	Gnb2	0.867	Uqcrc	0.574	CPR1	0.488	Ras	0.376	CS	0.840
11	PUP3	0.818	Cox5a	0.477	PUP3	0.520	Vps26	0.207	ATP2	0.703
12	RPT6	0.799	CPR1	0.643	Copb1	0.612	ATPDH	0.411	PUP3	0.785
13	NDUFB9	0.793	PUP3	0.592	CS	0.649	ATP2	0.249	ATP1	0.856
14	ATPDH	0.966	RPT6	0.643	ATP1	0.647	NdufS7	0.231	Ras	0.863
15	CPR1	0.804	CS	0.542	Uqcr8	0.589	Uqcr8	0.180	Uqcr8	0.733
16	PsmB6	0.894	ATP2	0.462	Gnb2	0.711	NDUFB9	0.244	α-TUB	0.777
17	SPRYp	0.763	ATP1	0.504	PsmB6	0.518	NdufA2	0.212	ATPDH	0.822
18	β-TUB2	0.766	Gnb2	0.536	ATP2	0.470	FKS1	0.277	CPR1	0.670
19	NdufA2	0.380	SPRYp	0.636	Cox5a	0.707	ATP1	0.340	Gnb2	0.823
20	ATP7	1.009	FKS1	0.463	SPRYp	0.595	Copb1	0.176	NDUFB9	0.800
21	Copb1	1.256	PsmB6	0.647	Vps26	0.554	RPT6	0.261	RPT6	0.737
22	Uqcrc	1.152	Ras	0.554	RPT6	0.593	Gnb2	0.319	PsmB6	0.709
23	Cox5a	1.008	ATPDH	0.531	NDUFB9	0.701	PsmB6	0.347	SPRYp	0.697
24	ATP1	1.043	β-TUB2	0.563	GAPDH	0.660	PUP3	0.129	NdufA2	0.451
25	Ras	1.035	NdufA2	0.225	β-TUB2	0.996	ATP7	0.160	β-TUB2	1.059
26	GAPDH	1.149	GAPDH	1.089	NdufA2	0.445	GAPDH	1.094	GAPDH	1.026

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.

[figures omitted; refer to PDF]

3.4. 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 ($\text{Vn}/\text{Vn}+1$) in each sample set. As shown in Figure 4, all the $\text{Vn}/\text{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.

[figure omitted; refer to PDF]

3.5. 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 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.

[figures omitted; refer to PDF]

4. 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–48]. Most of them, however, investigate the stability of traditional housekeeping genes such as GAPDH, actin, 18S rRNA, tubulins, 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 NormFinder 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 ($\text{Vn}/\text{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 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.

Acknowledgments

This work was supported by Public Welfare Industry (Agriculture) Scientific Research Project (201503137), Fundamental Research Project for Public Welfare Scientific Research Institutes in Fujian (2017R1020-1), Seed Industry Innovation and Industrialization Project of Fujian Province (fjzycxny2017009), and Fundamental Research Project for Public Welfare Scientific Research Institutes in Fujian (2019R1035-2).

References

[1] R. Narsai, A. Ivanova, S. Ng, J. Whelan, "Defining reference genes in Oryza sativa using organ, development, biotic and abiotic transcriptome datasets," BMC Plant Biology, vol. 10 no. 1,DOI: 10.1186/1471-2229-10-56, 2010.

[2] S. Dixit, V. K. Jangid, A. Grover, "Evaluation of suitable reference genes in Brassica juncea and its wild relative Camelina sativa for qRT-PCR analysis under various stress conditions," PLoS One, vol. 14 no. 9, article e0222530,DOI: 10.1371/journal.pone.0222530, 2019.

[3] K. Niu, Y. Shi, H. Ma, "Selection of candidate reference genes for gene expression analysis in Kentucky Bluegrass ( Poa pratensis L.) under abiotic stress," Frontiers in plant science, vol. 8,DOI: 10.3389/fpls.2017.00193, 2017.

[4] T. Zheng, Z. Chen, Y. Ju, H. Zhang, M. Cai, H. Pan, Q. Zhang, "Reference gene selection for qRT-PCR analysis of flower development in Lagerstroemia indica and L. speciosa," PloS one, vol. 13 no. 3,DOI: 10.1371/journal.pone.0195004, 2018.

[5] W. Li, L. Zhang, Y. Zhang, G. Wang, D. Song, Y. Zhang, "Selection and validation of appropriate reference genes for quantitative real-time PCR normalization in staminate and perfect flowers of andromonoecious Taihangia rupestris," Frontiers in Plant Science, vol. 8,DOI: 10.3389/fpls.2017.00729, 2017.

[6] M. Chen, J. Liao, H. Li, Z. Cai, Z. Guo, M. P. Wach, Z. Wang, "iTRAQ-MS/MS proteomic analysis reveals differentially expressed proteins during post-harvest maturation of the white button mushroom Agaricus bisporus," Current Microbiology, vol. 74 no. 5, pp. 641-649, DOI: 10.1007/s00284-017-1225-y, 2017.

[7] F. Li, L. D. Chen, L. Y. Ai, Y. C. Liu, M. Yan, S. J. Sun, "Comparative transcriptomics analyses of Pleurotus eryngii at different developmental stages," Mycosystema, vol. 37 no. 12, pp. 1586-1597, 2018.

[8] M. N. Heneghan, C. Burns, A. M. Costa, K. S. Burton, M. P. Challen, A. M. Bailey, G. D. Foster, "Functional analysis of Agaricus bisporus serine proteinase 1 reveals roles in utilization of humic rich substrates and adaptation to the leaf-litter ecological niche," Environmental Microbiology, vol. 18 no. 12, pp. 4687-4696, DOI: 10.1111/1462-2920.13350, 2016.

[9] R. Hu, C. Fan, H. Li, Q. Zhang, Y. F. Fu, "Evaluation of putative reference genes for gene expression normalization in soybean by quantitative real-time RT-PCR," BMC Molecular Biology, vol. 10 no. 1,DOI: 10.1186/1471-2199-10-93, 2009.

[10] X. Niu, J. Qi, G. Zhang, J. Xu, A. Tao, P. Fang, J. Su, "Selection of reliable reference genes for quantitative real-time PCR gene expression analysis in Jute ( Corchorus capsularis ) under stress treatments," Frontiers in Plant Science, vol. 6,DOI: 10.3389/fpls.2015.00848, 2015.

[11] J. Li, J. Han, Y. Hu, J. Yang, "Selection of reference genes for quantitative real-time PCR during flower development in tree peony ( Paeonia suffruticosa Andr.)," Frontiers in plant science, vol. 7 no. 7,DOI: 10.3389/fpls.2016.00516, 2016.

[12] Y. Zhang, X. Han, S. Chen, L. Zheng, X. He, M. Liu, G. Qiao, Y. Wang, R. Zhuo, "Selection of suitable reference genes for quantitative real-time PCR gene expression analysis in Salix matsudana under different abiotic stresses," Scientific Reports, vol. 7 no. 1, 2017.

[13] L. Gutierrez, M. Mauriat, S. Guenin, J. Pelloux, J. Lefebvre, R. Louvet, C. Rusterucci, T. Moritz, F. Guerineau, C. Bellini, O. Van Wuytswinkel, "The lack of a systematic validation of reference genes: a serious pitfall undervalued in reverse transcription-polymerase chain reaction (RT-PCR) analysis in plants," Plant Biotechnology Journal, vol. 6 no. 6, pp. 609-618, DOI: 10.1111/j.1467-7652.2008.00346.x, 2008.

[14] A. Kundu, A. Patel, A. Pal, "Defining reference genes for qPCR normalization to study biotic and abiotic stress responses in Vigna mungo," Plant Cell Reports, vol. 32 no. 10, pp. 1647-1658, DOI: 10.1007/s00299-013-1478-2, 2013.

[15] H. Li, W. Yao, Y. Fu, S. Li, Q. Guo, "De novo assembly and discovery of genes that are involved in drought tolerance in Tibetan Sophora moorcroftiana," PLoS One, vol. 10 no. 1,DOI: 10.1371/journal.pone.0111054, 2015.

[16] H. Muraguchi, K. Umezawa, M. Niikura, M. Yoshida, T. Kozaki, K. Ishii, K. Sakai, M. Shimizu, K. Nakahori, Y. Sakamoto, C. Choi, C. Y. Ngan, E. Lindquist, A. Lipzen, A. Tritt, S. Haridas, K. Barry, I. V. Grigoriev, P. J. Pukkila, "Strand-specific RNA-seq analyses of fruiting body development in Coprinopsis cinerea," PLoS One, vol. 10 no. 10, article e0141586,DOI: 10.1371/journal.pone.0141586, 2015.

[17] Q. Zhang, J. Zhang, M. Gong, R. Pan, Y. Liu, L. Tao, K. He, "Transcriptome analysis of the gene expression profiles associated with fungal keratitis in mice based on RNA-Seq," Investigative Ophthalmology & Visual Science, vol. 61 no. 6, pp. 32-32, DOI: 10.1167/iovs.61.6.32, 2020.

[18] J. Zhang, A. Ren, H. Chen, M. Zhao, L. Shi, M. Chen, H. Wang, Z. Feng, "Transcriptome analysis and its application in identifying genes associated with fruiting body development in basidiomycete Hypsizygus marmoreus," PLoS One, vol. 10 no. 4, article e0123025,DOI: 10.1371/journal.pone.0123025, 2015.

[19] W. Zheng, T. Peng, W. He, H. Zhang, "High-throughput sequencing to reveal genes involved in reproduction and development in Bactrocera dorsalis (Diptera: Tephritidae)," PloS one, vol. 7 no. 5,DOI: 10.1371/journal.pone.0036463, 2012.

[20] H. Yang, J. Liu, S. Huang, T. Guo, L. Deng, W. Hua, "Selection and evaluation of novel reference genes for quantitative reverse transcription PCR (qRT-PCR) based on genome and transcriptome data in Brassica napus L," Gene, vol. 538 no. 1, pp. 113-122, DOI: 10.1016/j.gene.2013.12.057, 2014.

[21] H. Zhuang, Y. Fu, W. He, L. Wang, Y. Wei, "Selection of appropriate reference genes for quantitative real-time PCR in Oxytropis ochrocephala Bunge using transcriptome datasets under abiotic stress treatments," Frontiers in Plant Science, vol. 6, 2015.

[22] F. Firenzuoli, L. Gori, G. Lombardo, "The medicinal mushroom Agaricus blazei murrill: review of literature and pharmaco-toxicological problems," Evidence-based Complementary and Alternative Medicine, vol. 5 no. 1,DOI: 10.1093/ecam/nem007, 2008.

[23] H. Wang, Z. Fu, C. Han, "The medicinal values of culinary-medicinal royal sun mushroom ( Agaricus blazei Murrill)," Evidence-Based Complementary and Alternative Medicine, vol. 2013,DOI: 10.1155/2013/842619, 2013.

[24] L. Wang, H. Li, H. Wei, X. Wu, L. Ke, "Identification of cadmium-induced Agaricus blazei genes through suppression subtractive hybridization," Food and chemical toxicology, vol. 63, pp. 84-90, DOI: 10.1016/j.fct.2013.10.036, 2014.

[25] Y. P. Lu, J. H. Liao, Z. J. Guo, Z. X. Cai, M. Y. Chen, "Genome survey and transcriptome analysis on mycelia and primordia of Agaricus blazei," BioMed Research International, vol. 2020 no. 11,DOI: 10.1155/2020/5204348, 2020.

[26] Y. P. Lu, Z. J. Guo, Z. X. Cai, M. Y. Chen, J. H. Liao, "Transcriptome analysis of Agaricus blazei fruiting bodies at different developmental stages," Mycosystema, vol. 38 no. 12, pp. 2161-2173, 2019.

[27] H. J. M. de Jonge, R. S. N. Fehrmann, E. S. J. M. de Bont, R. M. W. Hofstra, F. Gerbens, W. A. Kamps, E. G. E. de Vries, A. G. J. van der Zee, G. J. te Meerman, A. ter Elst, "Evidence based selection of housekeeping genes," PLoS One, vol. 2 no. 9,DOI: 10.1371/journal.pone.0000898, 2007.

[28] T. Czechowski, M. Stitt, T. Altmann, M. K. Udvardi, W. R. Scheible, "Genome-wide identification and testing of superior reference genes for transcript normalization in Arabidopsis," Plant Physiology, vol. 139 no. 1,DOI: 10.1104/pp.105.063743, 2005.

[29] Y. Tao, A. F. Van Peer, Q. Huang, Y. Shao, L. Zhang, B. Xie, B. Xie, "Identification of novel and robust internal control genes from Volvariella volvacea that are suitable for RT-qPCR in filamentous fungi," Scientific Reports, vol. 6 no. 1,DOI: 10.1038/srep29236, 2016.

[30] S. Lalitha, "Primer premier 5," Biotech Software & Internet Report: The Computer Software Journal for Scient, vol. 1 no. 6, pp. 270-272, DOI: 10.1089/152791600459894, 2000.

[31] Q. Q. Tan, L. Zhu, Y. Li, W. Liu, W. H. Ma, C. L. Lei, X. P. Wang, "A de novo transcriptome and valid reference genes for quantitative real-time PCR in Colaphellus bowringi," PLoS One, vol. 10 no. 2, article e0118693,DOI: 10.1371/journal.pone.0118693, 2015.

[32] T. Arima, M. Yamamoto, A. Hirata, S. Kawano, T. Kamada, "The eln3 gene involved in fruiting body morphogenesis of Coprinus cinereus encodes a putative membrane protein with a general glycosyltransferase domain," Fungal Genetics and Biology, vol. 41 no. 8, pp. 805-812, DOI: 10.1016/j.fgb.2004.04.003, 2004.

[33] K. J. Livak, T. D. Schmittgen, "Analysis of relative gene expression data using real-time quantitative PCR and the 2 -△△CT method," Methods, vol. 25 no. 4, pp. 402-408, DOI: 10.1006/meth.2001.1262, 2001.

[34] J. Vandesompele, "Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes," Genome biology, vol. 3 no. 7, article research0034.1,DOI: 10.1186/gb-2002-3-7-research0034, 2002.

[35] C. L. Andersen, J. L. Jensen, T. F. Orntoft, "Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets," Cancer Research, vol. 64 no. 15, pp. 5245-5250, DOI: 10.1158/0008-5472.CAN-04-0496, 2004.

[36] J. Y. Wu, B. He, Y. J. Du, W. C. Li, Y. Z. Wei, "Analysis method of systematically evaluating stability of reference genes using geNorm, NormFinder and BestKeeper," Modern Agricultural Science and Technology, vol. 5, pp. 278-281, 2017.

[37] Z. Yang, Y. Chen, B. Hu, Z. Tan, B. Huang, "Identification and validation of reference genes for quantification of target gene expression with quantitative real-time PCR for tall fescue under four abiotic stresses," PLoS One, vol. 10 no. 3, 2015.

[38] M. W. Pfaf, A. Tichopad, C. Prgomet, T. P. Neuvians, "Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper – Excel-based tool using pair-wise correlations," Biotechnology Letters, vol. 26 no. 6, pp. 509-515, DOI: 10.1023/B:BILE.0000019559.84305.47, 2004.

[39] X. Zhu, X. Li, W. Chen, J. Chen, W. Lu, L. Chen, D. Fu, "Evaluation of new reference genes in papaya for accurate transcript normalization under different experimental conditions," PLoS ONE, vol. 7 no. 8, article e44405,DOI: 10.1371/journal.pone.0044405, 2012.

[40] E. Chang, S. Shi, J. Liu, T. Cheng, L. Xue, X. Yang, W. Yang, Q. Lan, Z. Jiang, "Selection of reference genes for quantitative gene expression studies in Platycladus orientalis (Cupressaceae) using real-time PCR," PLoS One, vol. 7 no. 3, article e33278,DOI: 10.1371/journal.pone.0033278, 2012.

[41] X. Xiao, J. Ma, J. Wang, X. Wu, P. Li, Y. Yao, "Validation of suitable reference genes for gene expression analysis in the halophyte Salicornia europaea by real-time quantitative PCR," Frontiers in Plant Science, vol. 5,DOI: 10.3389/fpls.2014.00788, 2015.

[42] R. Castanera, L. López-Varas, A. G. Pisabarro, L. Ramírez, "Validation of reference genes for transcriptional analyses in Pleurotus ostreatus by using reverse transcription-quantitative PCR," Applied and Environmental Microbiology, vol. 81 no. 12, pp. 4120-4129, DOI: 10.1128/AEM.00402-15, 2015.

[43] N. V. Demidenko, M. D. Logacheva, A. A. Penin, "Selection and validation of reference genes for quantitative real-time PCR in buckwheat ( Fagopyrum esculentum ) based on transcriptome sequence data," PLoS One, vol. 6 no. 5,DOI: 10.1371/journal.pone.0019434, 2011.

[44] S. Artico, S. M. Nardeli, O. Brilhante, M. F. Grossi-de-Sa, M. Alves-Ferreira, "Identification and evaluation of new reference genes in Gossypium hirsutum for accurate normalization of real-time quantitative RT-PCR data," BMC Plant Biology, vol. 10 no. 1, 2010.

[45] M. Delporte, G. Legrand, J. L. Hilbert, D. Gagneul, "Selection and validation of reference genes for quantitative real-time PCR analysis of gene expression in Cichorium intybus," Frontiers in Plant Science, vol. 6,DOI: 10.3389/fpls.2015.00651, 2015.

[46] Y. Chen, B. Hu, Z. Tan, J. Liu, Z. Yang, Z. Li, B. Huang, "Selection of reference genes for quantitative real-time PCR normalization in creeping bentgrass involved in four abiotic stresses," Plant Cell Reports, vol. 34 no. 10, pp. 1825-1834, DOI: 10.1007/s00299-015-1830-9, 2015.

[47] J. Qian, Y. Gao, Y. Wáng, Y. Wu, Y. Wāng, Y. Zhao, H. Chen, D. Bao, J. Xu, X. Bin, "Selection and evaluation of appropriate reference genes for RT-qPCR normalization of Volvariella volvacea gene expression under different conditions," BioMed research international, vol. 2018,DOI: 10.1155/2018/6125706, 2018.

[48] C. Yang, L. Ma, D. Xiao, Z. Ying, X. Jiang, Y. Lin, "Identification and evaluation of reference genes for qRT-PCR normalization in Sparassis latifolia (Agaricomycetes)," International journal of medicinal mushrooms, vol. 21 no. 3, 2019.

[49] F. Jacob, R. Guertler, S. Naim, S. Nixdorf, A. Fedier, N. F. Hacker, V. Heinzelmann-Schwarz, "Careful selection of reference genes is required for reliable performance of RT-qPCR in human normal and cancer cell lines," PLoS One, vol. 8 no. 3, article e59180,DOI: 10.1371/journal.pone.0059180, 2013.

[50] M. Kanakachari, A. Solanke, N. Prabhakaran, I. Ahmad, G. Dhandapani, N. Jayabalan, P. A. Kumar, "Evaluation of suitable reference genes for normalization of qPCR gene expression studies in brinjal ( Solanum melongena L.) during fruit developmental stages," Applied Biochemistry and Biotechnology, vol. 178 no. 3, pp. 433-450, DOI: 10.1007/s12010-015-1884-8, 2016.

[51] S. Qi, L. Yang, X. Wen, Y. Hong, X. Song, M. Zhang, S. Dai, "Reference gene selection for RT-qPCR analysis of flower development in Chrysanthemum morifolium and Chrysanthemum lavandulifolium," Frontiers in Plant Science, vol. 7,DOI: 10.3389/fpls.2016.00287, 2016.

[52] K. Dudziak, M. Sozoniuk, H. Szczerba, A. Kuzdraliński, K. Kowalczyk, A. Börner, M. Nowak, "Identification of stable reference genes for qPCR studies in common wheat ( Triticum aestivum L.) seedlings under short-term drought stress," Plant Methods, vol. 16,DOI: 10.1186/s13007-020-00601-9, 2020.

[53] X. Zhou, X. Wu, M. Chu, C. Liang, X. Ding, J. Pei, L. Xiong, P. Bao, X. Guo, P. Yan, "Validation of suitable reference genes for gene expression studies on yak testis development," Animals, vol. 10 no. 2,DOI: 10.3390/ani10020182, 2020.

[54] K. E. Reid, N. Olsson, J. Schlosser, F. Peng, S. T. Lund, "An optimized grapevine RNA isolation procedure and statistical determination of reference genes for real-time RT-PCR during berry development," BMC Plant Biology, vol. 6 no. 1,DOI: 10.1186/1471-2229-6-27, 2006.

[55] M. Expósito-Rodríguez, A. A. Borges, A. Borges-Pérez, J. A. Pérez, "Selection of internal control genes for quantitative real-time RT-PCR studies during tomato development process," BMC Plant Biology, vol. 8 no. 1,DOI: 10.1186/1471-2229-8-131, 2008.