BioMed Research International
Volume 2021 (2021), Article ID 6661842, 13 pages
https://doi.org/10.1155/2021/6661842

Research Article

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

Yuan-Ping Lu, Jian-Hua Liao, Zhong-Jie Guo, Hui-Qing Zheng, Ling-Fang Lu, Zhi-Xin Cai, Zhi-Heng Zeng, Zheng-He Ying, and Mei-Yuan Chen

Institute of Edible Fungi, Fujian Academy of Agricultural Sciences, Fuzhou, 350014 Fujian Province, China

Correspondence should be addressed to Jian-Hua Liao; ljhua1@sohu.com and Mei-Yuan Chen; cmy1972@gmail.com

Received 9 November 2020; Accepted 3 February 2021; Published 9 April 2021

Academic Editor: Giuliana Banche

Copyright © 2021 Yuan-Ping Lu et al. 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.

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.

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 ( and ) 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 () 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 (). 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 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 () [37]. We then imported relative quantities into geNorm and calculated the expression stability ( value) of each gene. value less than 1.5 indicated that the gene can be selected for a reference gene, and the gene with the lowest value was considered as the most stable gene [36, 37]. And then, pairwise variation () 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 [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 () 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-3)
Product size (bp)
Amplification efficiency (%)

Uqcr8
A08809
Ubiquinol-cytochrome c reductase complex subunit 8
F:CACCCAATATGCCATTTCTCC
R:AATCCCGTATCCGACACCAA
136
99.2
1.000

Vps26
A09130
Vacuolar protein sorting 26
F:GGGAGAAGATGCTCGGAAAC
R:CCGCCAACTCTTGCGATAA
230
90.0
0.998

Gnb2
A04205
Guanine nucleotide-binding protein beta subunit 2
F: ACAAAGTCCACGCCATCCCT
R: GAGACCACCGCAAGCAACAT
89
92.9
0.998

Copb1
A08816
Coatomer subunit beta
F:TTATTGGGAAGTCTGTCCGAAGT
R: TGAACGACAGGTCGGTATGAG
190
93.7
0.998

Uqcrc
A01805
Core subunit of the ubiquinol-cytochrome c reductase complex
F: ACAACCTCGCCAACTCATTCA
R:TGCTCATTCTTGCCCATTCTTT
141
101.4
0.997

ATP7
A05732
ATP7 subunit D of the stator stalk of mitochondrial F1F0 ATP synthase
F: TCGCTCATTATCGCACTATCCT
R:GGCACTGTCCATTTACCCTTCTT
319
93.7
0.999

NdufS7
A01315
NdufS7 ubiquinone oxidoreductase subunit
F: CCTTGAAGTCCGCAGTTTGTAATG
R: TTTGTTTGTGGAGACGAGGTGAAT
100
97.6
0.998

PUP3
A04280
20S proteasome subunit beta 3
F:GAAGAGCGTGAGATTGAACCAGA
R:TCTGTCGCAGCAATAAAGGGTT
152
95.2
0.996

PsmB6
A05770
Proteasome subunit beta type-6
F: CGTCCACGACCGCATTTACT
R:AGGTCCACCCGACATTTGC
112
96.8
0.999

Ras
A03682
Ras protein
F:GGCAACAAGTGCGACAAATCA
R:CCACAGCGAAAGCAACAGAAA
251
103.4
0.997

CPR1
A01938
Probable proteasome subunit beta type-4
F: CCTGCTGCTGCCTCGTCTT
R:CTCTTGAGTGTTTCTATGCCTTCC
260
92.7
0.996

CS
A02079
Citrate synthase
F:GCTTGTCGGTTCAGCCCTTTCT
R:TCTTCCTTAGAACGGCGTGGC
230
97.4
0.997

SPRYp
A04044
Endosomal SPRY domain protein
F: ACCCTATCCCAGTTCCTTTCCA
R:CTTCACCTCCCAGTAGTAGACTTCATT
198
103.0
0.998

GAPDH
A07487
Glyceraldehyde-3-phosphate dehydrogenase
F: TTCGGCTCATTTGAAGGGTG
R:TTGGTTGTGCAAGAAGCGTTAG
140
94.0
0.998

ATP2
A06017
ATP2 beta subunit of the F1 sector of mitochondrial F1F0 ATP synthase
F:CCGAATACTTCCGTGACGAGG
R:TTGGTAACCGACAGCAGAGGG
128
99.0
0.999

ATP1
A07930
ATP1 alpha subunit of the F1 sector of mitochondrial F1F0 ATP synthase
F: TCCCTCCTTCTCCGTCGTCC
R:TCTTAGTCTGGGCGGCTGAA
300
96.6
0.999

Cdc4
A08640
Myosin regulatory light chain cdc4
F: CAGGCAGAAGTCGCTGAGATTGT
R:GAGGATTGTGCGGACAAAGGAT
300
96.3
0.992

FKS1
A09827
1, 3-beta-glucan synthase
F: ACCATTCACGCCGATTACATT
R:GCATAGCATTACGCCAACGA
187
104.0
0.996

Cox5a
A03337
COX5A subunit VA of cytochrome c oxidase
F: TCTCGTCCGTTGTGCGTCTG
R:TCGCGGTAGCGTAGTCGTTGA
206
93.2
0.996

ATPDH
A06239
ATP-dependent helicase
F: TCAAATCTGGCGTTGTCCC
R:CCGATTCCTTGATGTTCATTACTC
211
104.3
0.986

α-TUB
A05905
Tubulin alpha
F:CATTTCCCACTTGCCACTTTC
R:CAACAGCAGCATTGACATCCTTA
202
100.3
0.999

NdufA2
A00266
NdufA2, NADH ubiquinone oxidoreductase 10.5kD subunit
F: GCGAACTCCGTTTCTTCTTCTCCC
R:AGTCTGGGTTGTGCTGTTTGATGG
102
96.2
0.994

RPT6
A01698
26S proteasome regulatory complex ATPase RPT6
F: GTTGTTTGTTATGGCTCGTGAA
R:AAGCAACTCCAACATCGTCCT
136
96.3
0.993

ATP3
A08690
ATP3 gamma subunit of the F1 sector of mitochondrial F1F0 ATP synthase
F: ACAAGTTCGTTTCTGCCATCTCC
R:TGGTAGCATCATCCTCGTTCTCAT
114
97.8
0.999

NDUFB9
A09944
NDUFB9 NADH-ubiquinone oxidoreductase complex I
F:TACCTCACAAACTCGCTTGACTGG
R:CGGTTACATCACGGTTCTTCTCG
100
92.8
0.994

β-TUB2
A00521
Tubulin beta-2
F:CGGGCAAGATGTAGAGCAC
R:TCCAGGATAAAAGGCGAAG
231
100.8
0.996

STT3
A07729
Oligosaccharyl transferase STT3 subunit
F: TGTAGCCAACACTGCCTCCT
R:TGCCAAAAGCACACAAATG
331
98.5
0.997

F and R indicated forward primer and reverse primer, respectively; 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 () of each primer pair ranged from 90.0% for Vps26 to 104.3% for ATPDH (Table 2). The correlation coefficient () 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 (RPT6) and minimum of (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 1: Cq values of 27 candidate reference genes in total samples of A. blazei. Squares indicated the mean Cq values in samples, and the bars indicated minimum to maximum values.

(a)

(b)

(c)

(d)

(e)

Figure 2: Expression stability of 26 candidate genes in all sample sets. Average expression stability () was calculated using geNorm. The least and the most stable genes were on the left and on the right, respectively.

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

values, negatively correlated with gene stability, were calculated using geNorm software. GeNorm analysis showed that all reference genes possessed an 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 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 (), and value. Genes with 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 and a significant value ( 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.

(a)

(b)

(c)

(d)

(e)

Figure 3: RefFinder-based ranking of the candidate reference genes. (a) Five developmental stages, (b) fruit body, (c) stipe, (d) pileus, and (e) total sample. The most stable genes were on the left and the least stable genes were on the right.

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 () in each sample set. As shown in Figure 4, all the 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 4: Determination of the optimal number of reference genes for normalization using geNorm. Pairwise variation () less than 0.15 indicated an additional reference gene was not required.

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.

(a)

(b)

(c)

(d)

Figure 5: Normalized expression of a target gene Ableln3. Relative normalization of Ableln3 in five developmental stages (a), fruit body (b), stipe (c), and pileus (d).

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 (), 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.

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.

Conflicts of Interest

The authors declare no conflict of interest.

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

Supplementary Materials

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, indicated DL2000 bp DNA marker. (b) The melting curves for 27 candidate reference genes and Ableln3 showed single peaks. (Supplementary Materials)

References

R. Narsai, A. Ivanova, S. Ng, and J. Whelan, “Defining reference genes in Oryza sativa using organ, development, biotic and abiotic transcriptome datasets,” BMC Plant Biology, vol. 10, no. 1, p. 56, 2010.

S. Dixit, V. K. Jangid, and 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, 2019.

K. Niu, Y. Shi, and 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, 2017.

T. Zheng, Z. Chen, Y. Ju et al., “Reference gene selection for qRT-PCR analysis of flower development in Lagerstroemia indica and L. speciosa,” PloS one, vol. 13, no. 3, p. e0195004, 2018.

W. Li, L. Zhang, Y. Zhang, G. Wang, D. Song, and 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, p. 729, 2017.

M. Chen, J. Liao, H. Li et al., “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, 2017.

F. Li, L. D. Chen, L. Y. Ai, Y. C. Liu, M. Yan, and S. J. Sun, “Comparative transcriptomics analyses of Pleurotus eryngii at different developmental stages,” Mycosystema, vol. 37, no. 12, pp. 1586–1597, 2018.

M. N. Heneghan, C. Burns, A. M. Costa et al., “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, 2016.

R. Hu, C. Fan, H. Li, Q. Zhang, and 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, p. 93, 2009.

X. Niu, J. Qi, G. Zhang et al., “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, 2015.

J. Li, J. Han, Y. Hu, and 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, 2016.

Y. Zhang, X. Han, S. Chen et al., “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, pp. 1–11, 2017.

L. Gutierrez, M. Mauriat, S. Guenin et al., “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, 2008.

A. Kundu, A. Patel, and 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, 2013.

H. Li, W. Yao, Y. Fu, S. Li, and 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, p. e111054, 2015.

H. Muraguchi, K. Umezawa, M. Niikura et al., “Strand-specific RNA-seq analyses of fruiting body development in Coprinopsis cinerea,” PLoS One, vol. 10, no. 10, article e0141586, 2015.

Q. Zhang, J. Zhang, M. Gong et al., “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, 2020.

J. Zhang, A. Ren, H. Chen et al., “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, 2015.

W. Zheng, T. Peng, W. He, and H. Zhang, “High-throughput sequencing to reveal genes involved in reproduction and development in Bactrocera dorsalis (Diptera: Tephritidae),” PloS one, vol. 7, no. 5, 2012.

H. Yang, J. Liu, S. Huang, T. Guo, L. Deng, and 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, 2014.

H. Zhuang, Y. Fu, W. He, L. Wang, and 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.

F. Firenzuoli, L. Gori, and 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, p. 15, 2008.

H. Wang, Z. Fu, and C. Han, “The medicinal values of culinary-medicinal royal sun mushroom (Agaricus blazei Murrill),” Evidence-Based Complementary and Alternative Medicine, vol. 2013, Article ID 842619, 2013.

L. Wang, H. Li, H. Wei, X. Wu, and L. Ke, “Identification of cadmium-induced Agaricus blazei genes through suppression subtractive hybridization,” Food and chemical toxicology, vol. 63, pp. 84–90, 2014.

Y. P. Lu, J. H. Liao, Z. J. Guo, Z. X. Cai, and M. Y. Chen, “Genome survey and transcriptome analysis on mycelia and primordia of Agaricus blazei,” BioMed Research International, vol. 2020, no. 11, p. 8, 2020.

Y. P. Lu, Z. J. Guo, Z. X. Cai, M. Y. Chen, and J. H. Liao, “Transcriptome analysis of Agaricus blazei fruiting bodies at different developmental stages,” Mycosystema, vol. 38, no. 12, pp. 2161–2173, 2019.

H. J. M. de Jonge, R. S. N. Fehrmann, E. S. J. M. de Bont et al., “Evidence based selection of housekeeping genes,” PLoS One, vol. 2, no. 9, 2007.

T. Czechowski, M. Stitt, T. Altmann, M. K. Udvardi, and W. R. Scheible, “Genome-wide identification and testing of superior reference genes for transcript normalization in Arabidopsis,” Plant Physiology, vol. 139, no. 1, pp. 5–17, 2005.

Y. Tao, A. F. Van Peer, Q. Huang et al., “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, 2016.

S. Lalitha, “Primer premier 5,” Biotech Software & Internet Report: The Computer Software Journal for Scient, vol. 1, no. 6, pp. 270–272, 2000.

Q. Q. Tan, L. Zhu, Y. Li et al., “A de novo transcriptome and valid reference genes for quantitative real-time PCR in Colaphellus bowringi,” PLoS One, vol. 10, no. 2, article e0118693, 2015.

T. Arima, M. Yamamoto, A. Hirata, S. Kawano, and 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, 2004.

K. J. Livak and 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, 2001.

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, 2002.

C. L. Andersen, J. L. Jensen, and 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, 2004.

J. Y. Wu, B. He, Y. J. Du, W. C. Li, and 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.

Z. Yang, Y. Chen, B. Hu, Z. Tan, and 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.

M. W. Pfaf, A. Tichopad, C. Prgomet, and 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, 2004.

X. Zhu, X. Li, W. Chen et al., “Evaluation of new reference genes in papaya for accurate transcript normalization under different experimental conditions,” PLoS ONE, vol. 7, no. 8, article e44405, 2012.

E. Chang, S. Shi, J. Liu et al., “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, 2012.

X. Xiao, J. Ma, J. Wang, X. Wu, P. Li, and 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, p. 788, 2015.

R. Castanera, L. López-Varas, A. G. Pisabarro, and 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, 2015.

N. V. Demidenko, M. D. Logacheva, and 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, 2011.

S. Artico, S. M. Nardeli, O. Brilhante, M. F. Grossi-de-Sa, and 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, pp. 1–12, 2010.

M. Delporte, G. Legrand, J. L. Hilbert, and 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, p. 651, 2015.

Y. Chen, B. Hu, Z. Tan et al., “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, 2015.

J. Qian, Y. Gao, Y. Wáng et al., “Selection and evaluation of appropriate reference genes for RT-qPCR normalization of Volvariella volvacea gene expression under different conditions,” BioMed research international, vol. 2018, 12 pages, 2018.

C. Yang, L. Ma, D. Xiao, Z. Ying, X. Jiang, and 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.

F. Jacob, R. Guertler, S. Naim et al., “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, 2013.

M. Kanakachari, A. Solanke, N. Prabhakaran et al., “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, 2016.

S. Qi, L. Yang, X. Wen et al., “Reference gene selection for RT-qPCR analysis of flower development in Chrysanthemum morifolium and Chrysanthemum lavandulifolium,” Frontiers in Plant Science, vol. 7, p. 287, 2016.

K. Dudziak, M. Sozoniuk, H. Szczerba et al., “Identification of stable reference genes for qPCR studies in common wheat (Triticum aestivum L.) seedlings under short-term drought stress,” Plant Methods, vol. 16, pp. 1–8, 2020.

X. Zhou, X. Wu, M. Chu et al., “Validation of suitable reference genes for gene expression studies on yak testis development,” Animals, vol. 10, no. 2, p. 182, 2020.

K. E. Reid, N. Olsson, J. Schlosser, F. Peng, and 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, p. 27, 2006.

M. Expósito-Rodríguez, A. A. Borges, A. Borges-Pérez, and 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, p. 131, 2008.