Validation of Internal Control Genes for Quantitative Real-Time PCR Gene Expression Analysis in Morchella

The reliability of qRT-PCR results depend on the stability of reference genes used for normalization, suggesting the necessity of identification of reference genes before gene expression analysis. Morels are edible mushrooms well-known across the world and highly prized by many culinary kitchens. Here, several candidate genes were selected and designed according to the Morchella importuna transcriptome data. The stability of the candidate genes was evaluated with geNorm and NormFinder under three different experimental conditions, and several genes with excellent stability were selected. The extensive adaptability of the selected genes was tested in ten Morchella species. Results from the three experimental conditions revealed that ACT1 and INTF7 were the most prominent genes in Morchella, CYC3 was the most stable gene in different development stages, INTF4/AEF3 were the top-ranked genes across carbon sources, while INTF3/CYC3 pair showed the robust stability for temperature stress treatment. We suggest using ACT1, AEF3, CYC3, INTF3, INTF4 and INTF7 as reference genes for gene expression analysis studies for any of the 10 Morchella strains tested in this study. The stability and practicality of the gene, vacuolar protein sorting (INTF3), vacuolar ATP synthase (INTF4) and14-3-3 protein (INTF7) involving the basic biological processes were validated for the first time as the candidate reference genes for quantitative PCR. Furthermore, the stability of the reference genes was found to vary under the three different experimental conditions, indicating the importance of identifying specific reference genes for particular conditions.


Primer Design and Verification
Based on the transcriptome sequencing results of M. importuna at different sclerotial developmental stages in our previous work, the reference genes expression primers were designed in this paper. Firstly, ten traditional reference genes of other crops were selected as candidate genes [39][40][41][42], such as lyceraldehyde-3-phosphate dehydrogenase (GPD), γ-actin (ACT), polyubiquitin (PUB), peptidyl-prolyl cis-trans isomerase (CYC), α-tubulin (ATU), elongation factor 1-α (AEF), HSP90 domain-containing protein (HSP90), Hsp70 protein-like protein (HSP70), and β-tubulin (BTU). Meanwhile, the relative expression of each gene at different developmental stages and the gene variation coefficient were obtained by quantitative analysis of the transcriptome data (All the transcriptome reads and final assembly results have been deposited at NCBI, under accession numbers SRR6655864-69 and GGGK00000000, respectively). The genes with the lowest variation coefficient were selected as candidate genes for further analysis [43]. A total of thirty-two reference genes were chosen and detected in this study.
Specific primers were designed using the NCBI Primer-BLAST tool. The information about gene names and primer sequence is provided in Table 1. The amplification products were assayed by electrophoresis in a 1% agarose gel, and those with a specific single band at an expected length were selected, thus thirteen interested genes were excluded due to minor specificity (data not shown).
The mean amplification efficiency of the remaining genes was estimated by relative quantitative PCR with a standard curve of four serial 10-fold dilution of cDNA (ranging from 100 ng to 100 pg) for measuring the Cq value to a specific threshold [44]. The amplification efficiency was calculated by the formula: E = 10ˆ(1/slope) [45]. For all the genes, the amplification efficiency values ranged from 93.3% to 107.9%, and the linear fit R 2 was greater than 0.98, except CYC3 with a high value of 117.4%. Simultaneously, one single peak was required in real-time melting curve analysis of each gene. Finally, nineteen genes were selected as the candidate reference genes for further evaluation (Table 1).

Verification and Expression Stability of Candidate Reference Genes
Expression levels of the candidate genes for all the samples in different development stages, carbon sources, and temperature stress are presented in Figure 1. The Cq values of candidate genes varied obviously in the three different experimental conditions, ranging from 15.81 of AEF3 to 32.25 of INTF1 ( Figure 1a). The stability of all the candidate genes was estimated by two Microsoft Excel-based tools. The nineteen genes in the fifteen samples were ranked using geNorm according to their expression stability measure M, based on the pairwise variation of a specific gene in all the reference genes. Consequently, the gene with the lowest M value was the most stable in expression [43], while the gene with the highest M value was poor in stability. The optimal number of reference genes for an accurate normalization was determined by the pairwise variation values (Vn/Vn+1) obtained by geNorm.

Verification and Expression Stability of Candidate Reference Genes
Expression levels of the candidate genes for all the samples in different development stages, carbon sources, and temperature stress are presented in Figure 1. The Cq values of candidate genes varied obviously in the three different experimental conditions, ranging from 15.81 of AEF3 to 32.25 of INTF1 (Figure 1a). The stability of all the candidate genes was estimated by two Microsoft Excel-based tools. The nineteen genes in the fifteen samples were ranked using geNorm according to their expression stability measure M, based on the pairwise variation of a specific gene in all the reference genes. Consequently, the gene with the lowest M value was the most stable in expression [43], while the gene with the highest M value was poor in stability. The optimal number of reference genes for an accurate normalization was determined by the pairwise variation values (Vn/Vn+1) obtained by geNorm. A cut-off value of 0.15 was used by Vandesompele et al. below which the inclusion of an additional reference gene is not required [43]. Compared with geNorm, NormFinder adds intra and inter group variations to calculation of normalization factors, and different results are expected due to distinct statistical algorithms of geNorm and NormFinder [43,46].
When considering all the datasets, ACT1/ATU1 (M = 0.54) was the best pair as internal control genes, while GPD2 was the least stable gene (M = 1.29) (Figure1b). However, INTF7/INTF4 was A cut-off value of 0.15 was used by Vandesompele et al. below which the inclusion of an additional reference gene is not required [43]. Compared with geNorm, NormFinder adds intra and inter group variations to calculation of normalization factors, and different results are expected due to distinct statistical algorithms of geNorm and NormFinder [43,46].
When considering all the datasets, ACT1/ATU1 (M = 0.54) was the best pair as internal control genes, while GPD2 was the least stable gene (M = 1.29) (Figure 1b). However, INTF7/INTF4 was determined by NormFinder as the most ideal pair of reference genes. In both programs, HP701, BTU1 and GPD2 were ranked as the least stable genes (Tables 2 and 3; Figure 1b).   (Tables 2 and 3; Figure 1c), indicating that they were not fit for quantitative real-time PCR normalization in M. importuna under such conditions.
For different carbon sources, INTF2 and AEF3 had an equivalent M value of 0.23, and were defined as the most stable pair of genes by geNorm. However, INTF4 was identified as an ideal gene by NormFinder, followed by AEF3, while INTF2 was ranked as the fourth reliable gene (Tables 2 and 3; Figure 1d). Meanwhile, CYC3 and ATU2 were excluded due to the lowest expression level in different carbon sources.
For temperature stress, CYC3/INTF3 showed the lowest M value of 0.19, followed by ATU1 (M = 0.23) whereas the ATU1/INTF3 pair showed the best performance as determined by NormFinder. In both geNorm and NormFinder software, INTF1 and HP701 showed the highest variation in all the internal control genes (Tables 2 and 3; Figure 1e).
Moreover, the pairwise variation values obtained by geNorm for the least number of reference genes in an accurate normalization were also assessed. All the V values for each experimental set were below 0.15, indicating that there is no need to add a third gene as a reference gene. However, considering the total datasets, the optimum number of reference genes increased to three (Figure 1f).

Validation of the Selected Reference Genes in Morchella
The multigene phylogenetic analysis showed that there are about 68 phylogenetic species in Morchella all over the world, which belong to three branches: the basic branch Rufubrunnea clade, the sister branches Esculenta clade and Elata clade [14,47,48]. The basic branch Rufubrunnea clade contains two species, Morchella anatolica and Morchella rufobrunnea, which occurred only in North America and Europe [47]. The Esculenta and Elata clade formed a large group of 66 different species. Recent studies have shown that there are 30 Morchella phylogenetic species attributed to the Esculenta and Elata clades in China [49]. In order to validate the stability and generality of the selected reference genes in different development stages in Morchella, the CYC3 and INF5 primer pair was selected and evaluated. The amplification results revealed that CYC3 showed stable expression, but INTF5 was very unstable in all the ten morel species (Figure 2d).  The top-ranked reference genes for different carbon sources and temperature treatments were also validated. INTF4 and AEF3 performed well in different carbon sources while INTF2 only showed ideal stability among eight species (Figure 2c,g,e). Additionally, INTF3 and ATU1 were two top-ranked genes in temperature treatment (Figure 2f,b). INTF3 showed high consistency in all the species and ATU1 performed well except for two species in Esculenta (Mes-6 and Mes-24) ( Table 4). Moreover, ACT1 and INTF7 were further examined due to their high stability in all the data (Figure 2a,h). In addition, the dissolution curves were used to confirm the amplification efficiency of primers in ten samples. The results showed that the selected primers (ACT1, AEF3, INTF4 and INTF7) had a single dissolution curve, and the results were consistent with the electrophoresis results (Figure 2a,c,g,h; Figure 3). The results confirmed that both of them could serve as suitable reference genes for normalization in the genus of Morchella.

qRT-PCR of Amylase Gene under Different Carbon Sources with Different Internal Control Genes
To verify the importance of accurate normalization, RNA was prepared from five different carbon sources (G: glucose, SD: sawdust, WG: wheat grain, RS: rice, MIX: containing equal proportions of rice straw, sawdust and wheat grain as in materials and experimental treatments mentioned) and subjected to qRT-PCR for analysis of one hypothetical starch degradation gene. Data obtained was then normalized with both stable (INTF4/AEF3) and unstable (GPD1/CYC3) internal control genes, as suggested by the aforementioned results. As M. importuna has a good ability to degrade starch [14], the amylase genes, Amy_G (primer: GATTGGGAGCTTTGAGCAAG/ GTGTTCTCTTCTGGGTTGGA, amplification product size: 185 bp, Tm: 60 °C, annotation as

qRT-PCR of Amylase Gene under Different Carbon Sources with Different Internal Control Genes
To verify the importance of accurate normalization, RNA was prepared from five different carbon sources (G: glucose, SD: sawdust, WG: wheat grain, RS: rice, MIX: containing equal proportions of rice straw, sawdust and wheat grain as in materials and experimental treatments mentioned) and subjected to qRT-PCR for analysis of one hypothetical starch degradation gene. Data obtained was then normalized with both stable (INTF4/AEF3) and unstable (GPD1/CYC3) internal control genes, as suggested by the aforementioned results. As M. importuna has a good ability to degrade starch [14], the amylase genes, Amy_G (primer: GATTGGGAGCTTTGAGCAAG/GTGTTCTCTTCTGGGTTGGA, amplification product size: 185 bp, Tm: 60 • C, annotation as amylo-1,6-glucosidase compared to the Nr XP_007781986 with E-value = 0), was selected to analyze the expression levels of different carbon sources. As shown in Figure 4, two stable internal control genes INTF4 and AEF3 were selected for control, and the expression patterns of Amy_G gene in the five matrices were consistent, which expression in WG and MIX was higher than that in RSand SD, and is the lowest in RS (downregulation to 24% contrast to G) and the highest in WG (1.35 times upregulation contrast to G). However, when the unstable GPD1 or CYC3 was used as a control, the pattern of the Amy_G gene expression was different, which had the lowest in RS (downregulation to 5.8%) and the highest in MIX (1.6 times upregulation) under the GPD1 as control, and had the lowest in RS (downregulation to 7.0%) and the highest in WG (0.57 times) under the CYC3 as control. Selecting different reference genes as internal control, the maximum differential multiple of the Amy_G gene expression was 4.26 times (RS under AEF3 and GPD1 as control). It is thus obvious that normalization with different internal control genes would result in completely different conclusion, as also reported by others [50,51]. genes INTF4 and AEF3 were selected for control, and the expression patterns of Amy_G gene in the five matrices were consistent, which expression in WG and MIX was higher than that in RSand SD, and is the lowest in RS (downregulation to 24% contrast to G) and the highest in WG (1.35 times upregulation contrast to G). However, when the unstable GPD1 or CYC3 was used as a control, the pattern of the Amy_G gene expression was different, which had the lowest in RS (downregulation to 5.8%) and the highest in MIX (1.6 times upregulation) under the GPD1 as control, and had the lowest in RS (downregulation to 7.0%) and the highest in WG (0.57 times) under the CYC3 as control.
Selecting different reference genes as internal control, the maximum differential multiple of the Amy_G gene expression was 4.26 times (RS under AEF3 and GPD1 as control). It is thus obvious that normalization with different internal control genes would result in completely different conclusion, as also reported by others [50,51].

Discussion
In order to exclude variations caused by RNA quantity, integrity, reverse transcription efficiency and other factors [39], it is imperative to correct the expression of reference genes during qRT-PCR analysis. However, there is no universal reference gene whose expression level remains constant in all conditions, thus finding a particular reference gene for qRT-PCR in certain conditions is significant [53]. Traditional reference genes such as Actin, β-tublin, GAPDH, EF1α and 18S rRNA, which were broadly used as stable internal control genes in human and other model creatures in PCR reactions, were found not always stable in different organs or conditions [28,32,[54][55][56]. The normalization of reference genes is really essential before quantitative real-time PCR in specific species or new research materials.

Discussion
In order to exclude variations caused by RNA quantity, integrity, reverse transcription efficiency and other factors [39], it is imperative to correct the expression of reference genes during qRT-PCR analysis. However, there is no universal reference gene whose expression level remains constant in all conditions, thus finding a particular reference gene for qRT-PCR in certain conditions is significant [53]. Traditional reference genes such as Actin, β-tublin, GAPDH, EF1α and 18S rRNA, which were broadly used as stable internal control genes in human and other model creatures in PCR reactions, were found not always stable in different organs or conditions [28,32,[54][55][56]. The normalization of reference genes is really essential before quantitative real-time PCR in specific species or new research materials.
Identification and evaluation of reference genes have been performed in many species [28,32,34,35,40,41,54,[57][58][59][60], but no such studies have been performed in Morchella, even in pezizales. In our previous work, the transcriptome of Morchella importuna No-4 strain was sequenced (unpublished), with a purpose to obtain several housekeeping genes in Morchella under different development stages, carbon sources and temperature stresses for qRT-PCR analysis. In this study, the candidate genes were evaluated using the two frequently used tools of geNorm [43] and NormFinder [46]. The expression stability of the nineteen target genes showed high variation in the three experimental sets. The performance of the prominent reference genes was found extremely inconsistent under different conditions, whereas the least stable genes could be circumvented by both programs. In this case, traditional reference genes such as GPD2, BTU1 and ATU were all excluded as unsuitable internal control genes. Additionally, HP701 and INTF1 were confirmed as variable genes under the temperature stress set. HP701 gene was annotated as Hsp70 protein-like protein, which was an important part of the cell's machinery for protein folding and favorable to protect cells from stress [61]. INTF1 was annotated as similar to DNA damage response protein wss1 [62] ( Table 1). Functionally, the expression stability of the HSP701 and INTF1 genes is more susceptible to temperature stress. Thereby, GPD2, BTU1, ATU2, INTF1 and HP701 were considered inappropriate for normalization of genes in Morchella under such a condition.
ACT1 and ATU1 were confirmed as the most suitable reference genes by geNorm among all the datasets, although they were not top-ranked by NormFinder (INTF7 was ranked as the best one, which was annotated as 14-3-3 protein, a family of conserved regulatory molecules and expressed in all eukaryotic cells [63]) in M. importuna. Actin (ACT1) was used as a universal reference gene in many crops. 14-3-3 protein (INTF7) has the ability to bind a multitude of functionally diverse signaling proteins, including kinases, phosphatases, and transmembrane receptors. Functionally, the INTF7 gene can be used as a stable reference gene in crops. α-Tubulin (ATU) was used as a universal reference gene in many fungi [64,65], but ATU1 was not the top preferred gene here. Based on the integrated data here and from other morel species, ACT1 and INTF7 were chosen as the prominent pair of genes for normalization in Morchella.
In different development stages, INTF5 (annotated as similar in function to ubiquitin carboxyl-terminal hydrolase 4) [66] and CYC3 were assumed to be the best pair of reference genes. Consistently, CYC3 revealed stable expression in all the ten morel species (Figure 2d). However, IINTF5 was found extremely reliable as an internal control gene in M. importuna but not in the other morel strains, probably due to the matching problem of primers.
For different carbon sources, AEF3/INTF2 and INTF4 primer pairs performed well as determined by both geNorm and NormFinder software. These three genes were also evaluated with all the ten morel strains. INTF2 was found unstable for Mes-19 and Mes-24. Elongation factor 1-alpha (AEF3) was commonly used as a universal reference gene in other crops [41,59]. Vacuolar ATP synthase (INTF4, annotated as vacuolar ATP synthase 16kDa proteolipid subunit) is required for an array of molecular processes, including protein sorting, zymogen activation, receptor-mediated endocytosis and synaptic vesicle proton gradient generation [67], but so far, no research has been performed on the reference gene for gene expression analysis. Accordingly, AEF3/INTF4 was determined as the best pair of reference genes for different carbon sources.
Under temperature stress, four top-ranked genes have similar stability by both geNorm and NormFinder software, but with a slight difference in position. CYC3/INTF3 was the best pair of reference genes by geNorm, while ATU1/INTF3 was the most ideal pair of reference genes by NormFinder. Besides, INTF5 was ranked the third by NormFinder and the fourth place by geNorm under temperature treatment. Therefore, INTF3 was determined by both two programs as a stable reference gene for normalization of genes under different temperature stresses. INTF3, vacuolar protein sorting/targeting protein 10, functions as a sorting receptor in the golgi compartment required for the intracellular sorting and delivery of soluble vacuolar proteins, which has reason to be a stable candidate gene based on its biological function of protein [68]. CYC3 appeared to be the most controversial gene with high stability in temperature stress and development stage tests, but exhibited the opposite results in the carbon source tests (Tables 2 and 3). Therefore, CYC3 was suggested as an internal control gene in temperature stress and development stage tests in Morchella.
The optimal number of reference genes varied slightly between different experimental conditions and total database, and two reference genes were required for normalization of gene expression analysis in this paper.
In summary, the stability of reference genes was validated under different experimental conditions, demonstrating the importance of identification of specific reference genes prior to quantitative real-time PCR. Collectively, ACT1/INTF7, INTF4/AEF3, and INTF3/CYC3 were ranked as the most stable reference genes in Morchella based on the results from all the three experimental sets of different development stages, carbon sources and temperature stresses. This study has provided several ideal reference genes for further evaluation of quantitative real-time PCR gene expression analysis in Morchella.

Materials and Experimental Treatments
The commercial cultivation strain Morchella importuna No-4 was used in this study, which has been cultured in China for many years and stored in the Institute of Applied Mycology of Huazhong Agricultural University. The transcriptome of this strain at different sclerotial development stages has been sequenced and analyzed (data not shown), which allowed us to obtain a large number of candidate genes for the expression stability analysis. Fifteen samples of M. importuna were collected and divided into three experimental sets, including six development stages, six different carbon sources for mycelial growth and three temperature treatment sets.
To detect the gene expression during different development stages, six samples were obtained as mycelium, white sclerotia, mature sclerotia, primordium, small ascocarps and mature ascocarps. Furthermore, the mycelium of M. importuna was inoculated on six different carbon sources, including glucose, sawdust, wheat, rice, mix and a carbon-free CYM liquid medium, which were all with 2% final concentration improved on the basis of CYM liquid medium (yeast extracts 2 g L −1 , peptone 2 g L −1 , K 2 HPO 4 1 g L −1 , MgSO 4 0.5 g L −1 , KH 2 PO 4 0.46 g L −1 ). The glucose in CYM was separately replaced by sawdust, wheat and rice. The mix sample was a mixture of sawdust, wheat and rice with 0.67% final concentration for mycelium culture. The carbon source experiments were performed by incubation at 23 • C for seven days in the dark. The low and high temperature treatment sets were used as a temperature pressure test. The mycelia were incubated at 23 • C in CYM liquid medium for six days, then treated separately at 4 • C and 32 • C for 24 h, with untreated mycelium as a control. The as-treated mycelia were collected and rinsed in ddH 2 O, followed by dehydration through clean filter paper extrusion, freezing in liquid nitrogen and storage at −80 • C until RNA extraction.

Selection of Reference Genes and Primer Design
In our previous study, a transcriptome sequencing and de novo assembly work were conducted on M. importuna at different sclerotial development stages (under publication), thus enabling us to select and design the candidate reference genes in this paper. Meanwhile, the expression of the gene was analyzed (transcripts per million, TPM), and the variation coefficient was calculated for the transcript expression in the different stages [69]. According to the lowest coefficient of variation (CV < 0.08) and the size of expression abundance, the candidate genes were selected. The generally used housekeeping genes in other fungi, plants and animals were also employed based on the gene annotation results. The cDNA sequences of candidate genes were extracted, and the location information of introns and exons was determined through genome alignment from transcriptome to genome of M. importuna (genomic data has been submitted to the NCBI database, accession number: QOKS00000000). Primers were designed by using the NCBI primer-BLAST design tool (https://www.ncbi.nlm.nih.gov/tools/ primer-blast) [70]. The primer design principle was as follows: the size of the target fragment was between 150-300 bp, the temperature was at 60 ± 3 • C, the optimal length of primer was 20 bp, and the primers could cross one intron at least.
The specificity of all the designed primers was verified through 1% agarose gel electrophoresis before qRT-PCR validation. The amplification reaction of each pair of primers was performed with 0.3 µL of 4 × dNTPmix, 0.3 µL of rTaq (Takara, Nojihigashi 7-4-38, Kusatsu, Japan), 2 µL 10 × PCR buffer and 10 nM each of primer in a total volume of 20 µL, under the conditions of 95 • C for 3 min, followed by 35 cycles of 95 • C for 15 s, 60 • C for 15 s and 72 • C for 15 s. The size of the primer amplification product was consistent with the target fragment, with a single and clear band present under agarose gel electrophoresis, and such primers were selected for subsequent analysis.

Total RNA Isolation and cDNA Synthesis
The total RNA was extracted from 100 mg of each sample powder in liquid nitrogen, using the phenol/SDS method. DEPC-treated water was mixed with 10 × STE buffer and 10% sodium dodecyl sulfate at a volume ratio of 7/2/1 to release the nucleic acid, then the nucleic acid was extracted with phenol solution (phenol/chloroform/isoamyl alcohol at a volume ratio of 25/24/1). Isopropanol alcohol was used to deposit the total RNA, then the RNA deposit was dissolved in DEPC water, the residual DNA was digested by DNAase I (Thermo Scientific, Waltham, MA, USA) to remove DNA contamination, and the final pure RNA solution was stored at −80 • C. RNA concentration and quantity were determined using a DS-11+Spectrophotometer (Applied Denovix, Wilmington, DE, USA), and the RNA samples with only one single peak at 260 nm and the OD260/OD280 ratio between 1.8 and 2.2 were used for further analysis. The integrity of RNA was also assessed by electrophoresis on 1% agarose gel. The cDNA strand was synthesized using HiScript II Reverse Transcriptase (Vazyme, Nanjing, China) with 500 ng of total RNA in 20 µL reaction volume according to the protocol of the manufacturer. The cDNA products were diluted 1:5 with double-distilled water prior to quantitative real-time PCR analysis.

Quantitative Real-Time PCR Analysis
The amplification of qRT-PCR was carried out in a 96-well optical plate with a CFX Connect Real-Time System (Applied with Bio-Rad, Hercules, CA, USA). The reaction system consisted of 5 µL of SYBR Green Master Mix (AceQ TM qPCR, Vazyme), 2.5 nM of each primer, and 1 µL of diluted cDNA in a final volume of 10 µL. A two-step RT-PCR protocol was adopted for all amplifications: 95 • C for 5 min, followed by 42 cycles of 95 • C for 10 s and 60 • C for 30 s, the fluorescence was measured at the end of each cycle. After 42 cycles, a melting curve was acquired by heating the amplicon from 65 to 95 • C to confirm whether a single product was generated. Then, the threshold values (Cq) of each sample were estimated, and the data were exported into Excel format for further primer expression stability analysis. The specificity of the amplification was analyzed through melting curve analysis [71]. All assays were performed with three technical replicates.
To estimate the expression stability of reference genes, Cq values were converted into non-normalized relative quantities, and corrected by PCR efficiency using the formula Q = (E + 1) −∆Cq where E represents the amplification efficiency of each gene, and ∆Cq represents the Cq values minus the lowest Cq value of each gene [52]. These data were imported into two Microsoft Excel-based statistics tools, geNorm [43] and NormFinder [46], which were performed according to the manuals.

Validation of Reference Genes
To confirm the amplification efficiency and adaptability of selected reference genes in different Morchella species, ten different morel species including five Esculenta clade species and five Elata clade species were chosen to examine the adaptability of eight ideal reference genes according to the previous results. The strains were identified by multigene phylogenetic analysis (Table 4) [47], whose adaptability was validated by using the same reaction system and protocol for quantitative real-time PCR analysis. After amplification, the PCR products were examined by electrophoresis in a 1.8% agarose gel. One single and clear band in an expected length indicated that the sample can be used as a suitable species adaptability reference gene for qRT-PCR normalization in Morchella.
Author Contributions: Q.Z., W.L. and Y.B. conceived and designed the experiments. Q.Z., Y.C., A.-F.L. and W.L. prepared the experiment materials. W.L., and Q.Z. carried out the experiments and collected the data. Data analysis was carried out by Q.Z. and W.L. Q.Z. wrote the manuscript. W.L. and Y.B. provided intellectual input and revised the manuscript. All authors have read and approved the final manuscript.