A Comparative Transcriptome Analysis of Volvariella volvacea Identified the Candidate Genes Involved in Fast Growth at the Mycelial Growth Stage

The edible straw mushroom, Volvariella volvacea, is one of the most important cultivated mushrooms in tropical and sub-tropical regions. Strain improvement for V. volvacea is difficult because of the unknown mechanisms involved in its growth regulation and substrate utilization. A comparative physiological and transcriptomic study was conducted between two commercially available straw mushroom strains (v9 and v26) to explore their fast-growth regulation mechanism(s). The physiological study showed that V. volvacea v9 had a shorter growth cycle and higher biological efficiency (4% higher) than that in v26. At least 14,556 unigenes were obtained from the four cDNA libraries (two replicates per strain). Among them, the expression of 1597 unigenes was up-regulated while 1352 were down-regulated. Four heat-shock proteins were highly expressed in v9, showing that v9 has the better ability to handle stresses and/or environmental changes. Moreover, up to 14 putative transporter genes were expressed at a higher level in v9 than those in v26, implying that v9 has a better ability to transport nutrients or export xenobiotics efficiently. Our report allows to identify the candidate genes involved in the fast growth requirement of V. volvacea, which represents a valuable resource for strain improvement in this commercially important edible mushroom.


Introduction
Straw mushroom (Volvariella volvacea), also called Chinese mushroom, has a long history of cultivation in tropical and subtropical regions [1,2]. In China, at least 330,000 tons of straw mushroom were produced in 2010 (accounting for more than 80% of the global production). The yield ranked the fifth among the commercially cultivated mushrooms [3]. Besides its high nutrient value, V. volvacea also has medicinal importance, including antitumor polysaccharides, immunosuppressive proteins and immunomodulatory lectins [4,5]. Despite its high demand in the mushroom markets, the conversion

Comparations of the Growth Rates between V. volvacea Strain v9 and Strain v26
V. volvacea v9 had a bell-shaped fruit body with a gray-white color. It was noted that its basal surface contacted closely to the cultivated materials. However, the fruit body of strain v26 was oval-shaped with a gray color. The contact area between the fruiting body and the cultivated materials was smaller than that of v9 ( Figure 1A). On the PDA media, mycelial growth rates of v9 and v26 were 10.2 ± 2 mm/d and 8.3 ± 2 mm/d, respectively ( Figure 1B). The following formula for biological efficiency was used: biological efficiency (%) = fresh fruit body weight/substrate dry weight × 100. The biological efficiency of strain v9 was 26.3%, which was higher than strain v26 (22.3%) ( Figure 1C). Three development stage periods were measured (egg stage, elongation stage and mature stage) [18]. In the egg stage, the stipe was hidden and ovoid; in the elongation stage, the stipe was stretched out of the universal veil and the pileus was not opened; in the mature stage, the pileus was fully expanded. The egg stage, elongation stage and mature stage time in strain v9 was 21.9, 11.5 and 7.3 h, respectively. Instead, the mature time for strain v26 was 25.6, 13.6 and 10.5 h, respectively ( Figure 1D). Collectively, v9 had higher biological efficiency and a shorter mature time than v26.

Sequencing, Assembly and Functional Annotation of V. volvacea
The four cDNA libraries from the mycelial growth stage yielded a total of 67 million raw reads. Up to 16 million reads per sample were obtained after the data filtering and trimming (Table 1). A total of 14,556 unigenes were assembled for all samples (Table 1). Among them, at least 85% of the unigenes were annotated using the databases of NR and KEGG (Table 1). A total of 58% and 43% of the unigenes encoded the "hypothetical proteins" in strains v9 and v26, respectively. The above observation agreed that many genes in V. volvacea were not able to be assigned to the known functions in Basidiomycetes [7,9]. Both strains were cultivated on rice straw and cotton waste. (D) Developmental periods of v26 and v9. Strain v9 developed faster than v26. Each experiment was repeated three times. Mycelial growth rates, biological efficiency and developmental stages time were carried out by a Student's t test. The mycelial growth rates, biological efficiency and egg stage time were significantly different (p < 0.05, The elongation stage and mature stage developmental time were not significantly different (p > 0.05). Bar = 1 cm in A. Data are presented as the average ± SD in B, C and D.

Sequencing, Assembly and Functional Annotation of V. volvacea
The four cDNA libraries from the mycelial growth stage yielded a total of 67 million raw reads. Up to 16 million reads per sample were obtained after the data filtering and trimming (Table 1). A total of 14,556 unigenes were assembled for all samples (Table 1). Among them, at least 85% of the unigenes were annotated using the databases of NR and KEGG (Table 1). A total of 58% and 43% of the unigenes encoded the "hypothetical proteins" in strains v9 and v26, respectively. The above observation agreed that many genes in V. volvacea were not able to be assigned to the known functions in Basidiomycetes [7,9].
Most of the unigenes were enriched in "Metabolism" and "Genetic Information Processing" (Figure 2A, left panel) using the KEGG pathway enrichment analysis. Among the metabolism pathways, the amino sugar and nucleotide sugar metabolism were most abundant, followed by nitrogen metabolism, glutathione metabolism, lysine degradation and fructose and mannose metabolism (Figure 2A, right panel). Enrichment analysis of Gene Ontology (GO) showed that "Oxidation-reduction process" had the most abundant genes. Moreover, a large number of genes involved in the "Carbohydrate Metabolic Process" and "Proteolysis" were found in the "Biological . Strain v9 developed faster than v26. Each experiment was repeated three times. Mycelial growth rates, biological efficiency and developmental stages time were carried out by a Student's t test. The mycelial growth rates, biological efficiency and egg stage time were significantly different (p < 0.05, The elongation stage and mature stage developmental time were not significantly different (p > 0.05). Bar = 1 cm in (A). Data are presented as the average ± SD in (B-D). Most of the unigenes were enriched in "Metabolism" and "Genetic Information Processing" (Figure 2A, left panel) using the KEGG pathway enrichment analysis. Among the metabolism pathways, the amino sugar and nucleotide sugar metabolism were most abundant, followed by nitrogen Genes 2020, 11, 161 4 of 14 metabolism, glutathione metabolism, lysine degradation and fructose and mannose metabolism ( Figure 2A, right panel). Enrichment analysis of Gene Ontology (GO) showed that "Oxidation-reduction process" had the most abundant genes. Moreover, a large number of genes involved in the "Carbohydrate Metabolic Process" and "Proteolysis" were found in the "Biological Process" (Figure 2B), indicating that V. volvacea has a strong potential to utilize plant-derived culture substrates. substrates.

Gene Expression and Identification of DEGs
At least 13,488 (92.7%), 13,427 (92.2%), 13,200 (90.7%) and 12,955 (89%) transcripts were identified in v9-1, v9-2, v26-1 and v26-2 libraries, respectively. The cluster analysis showed that v9-1 was close to v9-2 and v26-1 was clustered with v26-2. Instead, v9-1 and v9-2 separated with v26-1 and v26-2. This observation indicated that the gene expression patterns were different in the two V. volvacea strains ( Figure 3A). qRT-PCR analysis of the selected genes was performed to validate the transcriptome data. The same gene expression patterns were consistent with those by the transcriptome analysis (Supplementary Table S1).  At least 2949 transcripts were shared between v9 and v26 at the vegetative growth stage. Among them, the expression of 1597 unigenes was up-regulated while 1352 were down-regulated ( Figure 3B). Up-regulated DEGs were assigned to Genetic Information Processing by the KEGG enrichment analysis of DGEs, followed by Environmental Information Processing, Cellular Processes, Amino Acid Metabolism and Carbohydrate Metabolism, as well as several others ( Figure 3C). Down-regulated DGEs were enriched in Genetic Information Processing, followed by Carbohydrate Metabolism, Amino Acid Metabolism, Enzyme Families, Environmental Information Processing and Cellular Processes ( Figure 3D).
A further significant enrichment analysis of the KEGG pathway showed that both up-regulated DGEs and down-regulated DGEs were significantly enriched on the pathways of Starch and Sucrose Metabolism, Thiamine Metabolism and ABC Transporters ( Figure 3E,F). Moreover, up-regulated DGEs were significantly enriched in the MAPK signaling pathway within Environmental Information Processing and Phagosome within cellular processes ( Figure 3E). However, down-regulated DGEs were significantly enriched in the Propanoate Metabolism within Carbohydrate Metabolism, and pathways related to Amino Acid Metabolism of Valine, Leucine and Isoleucine Degradation, Glycine, Serine and Threonine Metabolism as well as Cyanoamino Acid Metabolism ( Figure 3F).
Thiamine plays a vital role in the metabolism of glucose and energy production. At least six genes involved in "Thiamine Metabolism" pathways were investigated in detail because they were significantly enriched in the KEGG pathways ( Figure 3E,F). The expression level of thiamine synthase genes iscS (ctg4698_g1) was higher in v9 than in that in v26, demonstrating a more active vitamin B1 synthesis activity. The thiamine metabolism gene THI20 (ctg9033_g1) and THI4 (ctg9315_g1) in v9 showed a lower expression level than that in v26. The expression of other selected thiamine metabolism genes was only slightly different between the strains.
The carbohydrate utilization gene is expressed differently in the two strains (Table 2). Three genes encoding the cellulase 1,4-beta-cellobiosidase (CBH1 and CBH2) were upregulated in v9, and two of the three lytic polysaccharide monooxygenase genes (LPMOs) were expressed at a higher level in v9. Moreover, at least 8 β-glycosidase genes (bglX) were downregulated in v9 compared to v26. Two (ctg9698_g1 and ctg11581_g1) of the three selected α-amylase genes (AMY) showed a higher expression level in v9 than they did in v26 (Table 2).

String Analysis
A total of 262 highly up-regulated DEGs and 302 highly down-regulated DEGs were obtained using string analysis (FDR < 0.001). A total of 58 up-regulated unigenes was annotated to 48 and 55 down-regulated unigenes was annotated to 44 ( Figure 4). As shown in Figure 4A, up-regulated genes could be grouped into heat-shock proteins (HSPs), sugar transporters, ABC transporters, cytochrome P450, NAD(P)-binding domain and others ( Figure 4A). Four heat-shock proteins, including HSP26, HSP104, HSP42 and SIS1 (Type II HSP40 co-chaperone), were highly expressed in v9 (FPKM > 1400) while those in v26 were not (FPKM < 700). Fourteen putative transporters were expressed at a higher level in v9 than those in v26. These transporters seem to be involved in various nutrient transportation; for example, four sugar transporters were found: HXT17 (hexose transporter), SNF3 (low glucose sensor that regulates glucose transport), YFL040W (transporter member of the sugar porter family) and YHK8 (presumed antiporter of the major facilitator superfamily); four ABC transporters were found: ATM1 (exports mitochondrially synthesized precursors of iron-sulfur (Fe/S) clusters to the cytosol), PDR12 (weak-acid-inducible multidrug transporter), SNQ2 (multidrug transporter) and YBT1 (bile acid transport); as well as other transporters were found: BCS1 (AAA ATPase family), COX1 (subunit I of cytochrome c oxidase), RSB1 (transport sphingoid long chain base (LCB)), SMF2 (manganese transporter), TOM70 (transit peptide receptor) and TPO1 (polyamine transporter that recognizes spermine, putrescine, and spermidine). Moreover, six up-regulated and five down-regulated unigenes were significantly enriched in the pathway of the ABC transporters ( Table 4). The expression of gene ABCB1 (K05658, ATP-binding cassette, subfamily B (MDR/TAP), member 1) was higher in v9 than in v26. ABCB1 encoded the ATP-binding ATPases, which played important roles in the exportation of xenobiotics. Therefore, one can assume that v9 may more actively export xenobiotic molecules than v26.

String Analysis
A total of 262 highly up-regulated DEGs and 302 highly down-regulated DEGs were obtained using string analysis (FDR < 0.001). A total of 58 up-regulated unigenes was annotated to 48 and 55 down-regulated unigenes was annotated to 44 ( Figure 4). As shown in Figure 4A, up-regulated genes could be grouped into heat-shock proteins (HSPs), sugar transporters, ABC transporters, cytochrome P450, NAD(P)-binding domain and others ( Figure 4A). Four heat-shock proteins, including HSP26, HSP104, HSP42 and SIS1 (Type II HSP40 co-chaperone), were highly expressed in v9 (FPKM > 1400) while those in v26 were not (FPKM < 700). Fourteen putative transporters were expressed at a higher level in v9 than those in v26. These transporters seem to be involved in various nutrient transportation; for example, four sugar transporters were found: HXT17 (hexose transporter), SNF3 (low glucose sensor that regulates glucose transport), YFL040W (transporter member of the sugar porter family) and YHK8 (presumed antiporter of the major facilitator superfamily); four ABC transporters were found: ATM1 (exports mitochondrially synthesized precursors of iron-sulfur (Fe/S) clusters to the cytosol), PDR12 (weak-acid-inducible multidrug transporter), SNQ2 (multidrug transporter) and YBT1 (bile acid transport); as well as other transporters were found: BCS1 (AAA ATPase family), COX1 (subunit I of cytochrome c oxidase), RSB1 (transport sphingoid long chain base (LCB)), SMF2 (manganese transporter), TOM70 (transit peptide receptor) and TPO1 (polyamine transporter that recognizes spermine, putrescine, and spermidine). Moreover, six up-regulated and five down-regulated unigenes were significantly enriched in the pathway of the ABC transporters ( Table 4). The expression of gene ABCB1 (K05658, ATP-binding cassette, subfamily B (MDR/TAP), member 1) was higher in v9 than in v26. ABCB1 encoded the ATP-binding ATPases, which played important roles in the exportation of xenobiotics. Therefore, one can assume that v9 may more actively export xenobiotic molecules than v26. . STRING analysis. A protein interactive network was used to display the genes that were related to the highly up-regulated, differently expressed unigenes (A), or related to the highly downregulated, differently expressed unigenes (B). Gene products are represented with circles and known associations between each gene or gene product are represented with a connecting line. Nodes of genes participating in similar functions are circled in red or in blue or in green. The network cluster was based on the k-means clustering method. . STRING analysis. A protein interactive network was used to display the genes that were related to the highly up-regulated, differently expressed unigenes (A), or related to the highly down-regulated, differently expressed unigenes (B). Gene products are represented with circles and known associations between each gene or gene product are represented with a connecting line. Nodes of genes participating in similar functions are circled in red or in blue or in green. The network cluster was based on the k-means clustering method.

Discussion
The current cultivation mode of straw mushrooms are field cultivation and indoor cultivation [2,19,20]. Strain v9 is more appropriate for cultivation in the field while v26 is mainly cultured indoor for commercial production. Compared to the indoor cultivation, the filed cultivation mode encounters more dramatic temperature changes between day and night. Temperature fluctuations lead to a series of metabolism activity adaptions; for example, changes in unsaturated fatty acid biosynthesis, as well as the accumulation of trehalose and glycogen. The cold-shock and heat-shock proteins played important physiological roles during the temperature fluctuations periods [21,22]. It is well known that the production of heat shock proteins are promoted by abnormal stress factors in various fungi. Heat-shock proteins function as molecular chaperone that participate in refolding or degradation of stress-damaged proteins [23]. For example, exposure to low temperature (i.e., 4 • C) for more than eight hours will cause irreversible damage to the V. volvacea mycelium [24]. To respond to the dramatic temperature, a number of HSP genes in V. volvacea were down-regulated after low temperature exposure [1]. In this study, we also found that the expression levels of hsp104, hsp42 and hsp26 genes in v9 was much higher than those in v26 (Figure 4). High expression of heat-shock proteins in v9 may enhances its adaptability to temperature fluctuations. This may be one of the reasons why the v9 strain was selected for field cultivation The mycelium growth rate of the v9 strain is faster than that of the v26 strain, and the biological efficiency is higher (Figure 1). The mycelial growth performance and biological efficiency of a mushroom species is mainly attributed to its hydrolytic enzymes system [25,26]. V. volvacea had a higher number of enzymes related to the degradation of cellulose, hemicellulose and pectin compared with other basidiomycetes [3,20,27]. Our results showed several carbohydrate utilization genes in v9 was transcribed at higher levels than those in v26; for example, cellulase 1,4-β-cellobiosidase (CBH1 and CBH2), which may efficiently contribute to cellobiose releasing from cellulose. Cellobiose was an inducer for the expression of cellulase genes (e.g., 1,4-β-cellobiosidase genes) in many fungi, including T. reesei and V. volvacea [28,29]. LPMOs are involved in oxidatively cleaving glycosidic linkages, which leads the substrate to be more susceptible to hydrolysis by those conventional cellulases [30]. Further, β-glycosidase genes (bglX) were downregulated in v9. β-glycosidases release glucose from cellobiose, which represses the transcription of cellulase genes [31,32]. Such a regulation mechanism may enable v9 to maintain a higher level of expression of extracellular cellulase. Similarly, the β-D-xylosidase gene XYL4 was expressed at a higher level in v26 than it did in v9. β-D-xylosidases digested xylobiose to produce glucose, which is a repressor for xylanase gene expression [31,32]. Two (ctg9698_g1 and ctg11581_g1) of the three selected α-amylase genes (AMY) showed a higher expression level in v9 than they did in v26 ( Table 2), indicating that these two AMY gene was possibly involved in starch utilization in the fast-growth stage. The hexokinase gene HK (ctg4916_g1) transcription level in v9 was around 4-fold higher than in v26, indicating that v9 was able to metabolize and transform the simple soluble sugars (e.g., glucose or fructose) into energy for the fast growth. Moreover, two TREH genes for the conversion of trehalose to glucose were differently regulated (ctg9797_g3 was upregulated in v9). The transcriptional level of malZ showed a slightly higher expression in v26 than it did in v9, which indicated a faster conversion of maltose into glucose in v26 than in v9. The higher expression of these carbohydrate utilization genes may contribute to a faster mycelium growth rate and higher biological efficiency of v9.
Our quantitative analysis and comparison of gene expression profiles in two strains suggested that MAPK modules might be involved in mycelial growth in V. volvacea. However, it remains unclear what the exact role of these MAPK pathway genes are; identifying potential regulators and targets in MAPK pathways will provide insights into the signal transduction pathways of V. volvacea. Our data also revealed that gene expression profiles were different between v9 and v 26 when MAPK signaling pathway analysis was performed (Table 3). MAPK modules serve central roles in intracellular signal transduction and are evolutionary conserved in eukaryotic cells, including yeast and fungi [33,34]. The MAPK cascade are involved in mating and pathogenic development in Ustilago maydis [35]. In Lentinula edodes, the MAPK cascade was also reported to play important roles in light signal transmission and metabolic regulation [36].
We have found an array of genes that contribute to the growth of V. volvacea and should increase our understanding of the fundamental and cellular processes of V. volvacea for commercial production and industrial use.

Strains and Culture Conditions
The V. volvacea dikaryotic strain v26 was isolated in a greenhouse from Foshan Zhihua Fungus Co., Ltd. (Foshan, China), and v9 was isolated in a field from Zhongshan Minzhong Edible Fungus Cooperative; both were preserved in the Fungal Culture Collection of Guangdong Key Laboratory for New Technology Research of Vegetables. The microorganisms used in the current study were maintained on potato dextrose agar (PDA, Merck, Germany) at 30 • C. The mycelial growth rate was determined in a Petri dish with a diameter of 9 cm. A 0.5 cm 2 mycelial agar plug was transferred into the center of a fresh PDA plate and the colonies' radii were measured every 24 h. Straw mushroom strains were cultivated on mixed substrates including 50% rice straw and 50% waste cotton (each with 65% water content) at 30 • C. Various developmental stages of V. volvacea were differentiated using the methods described previously by Chang and Yau [18]. At least 10 replicates for each stage were analyzed.

RNA Extraction and cDNA Library Construction and Sequencing
V. volvacea strains v9 and v26 were cultured in 250 mL flasks containing 100 mL Potato Dextrose Broth (PDB) at 30 • C for 5 days. Mycelium were collected and immediately stored in liquid nitrogen. The two replicates per strain were used to construct the comparative cDNA libraries. Total RNA was extracted using Trizol reagents. A total of 1 µg of RNA was proceeded with the Illumina TruSeq RNA Sample Preparation kit v2, including steps to purify and fragment the mRNA, carry out first-strand cDNA synthesis and second-strand cDNA synthesis, repair the ends, adenylate the 3' ends and ligate the adaptors. Following PCR amplification, the libraries were validated using an Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA 95051, USA), which indicated average cDNA fragments of~360 bp in length.
Four cDNA libraries were sequenced using an Illumina MiSeq sequencing platform (available at Agro-biological Gene Research Center of Guangdong Academy of Agricultural Sciences, Guangzhou, China) by using a 2× 75 bp MiSeq reagent kit v3 (Illumina, San Diego, CA 92122, USA). Raw paired-end reads were submitted to the Sequence Read Archive of the NCBI (accession number: PRJNA408191).
To examine the expression of the 14 putative thiamine metabolism-as well as starch and sucrose metabolism-related genes, qRT-PCR assays were performed. The primers for the 14 genes related to thiamine metabolism, starch and sucrose metabolism are listed in Supplementary Table S1. The PCR mix was prepared with 1 µL cDNA samples as templates in the presence of a SYBR Green PCR Master Mix and gene-specific primers. The reactions were performed in a BIO-RAD Cycler IQ Multi-Color Real Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The relative levels of the amplified mRNA were evaluated according to the 2 −∆∆Ct method using the actin gene for normalization.

Sequence Processing and Bioinformatics Analysis
To obtain high quality sequences, the QC Toolkit (v2.3.1) was adopted to remove adaptor sequences, ambiguous reads and low-quality reads [37]. The clean reads from all samples were pooled together to perform the de novo assembly using the Trinity program (r20140717), which included three modules, namely inchworm, chrysalis and butterfly [38], with parameters of "-seqType fq -JM 50G -CPU 12 -trimmomatic". Then the nucleotide sequences of the unigenes were aligned to the databases of NR and KEGG using the Blastx program with the e-value lower than 1e-5. Moreover, the clean reads per sample were aligned to the assembled transcripts with using Bowtie [39]. RSEM [40] was used to calculate the FPKM values (Fragments Per kb per Million reads) to obtain gene expression level, which were done through a perl script (align_and_estimate_abundance.pl -transcripts Trinity.fasta -left SeqR1.fq -right SeqR2.fq -seqType fq -est_method RSEM -aln_method bowtie -trinity_mode -prep_reference). EdgeR [41] was performed to detect differential expression with an FDR ≤ 0.05 and a relative change threshold of 2-fold. KOBAS [42] was used for enrichment analysis of the Gene Ontology and KEGG pathway. The functional differences between the two strains were considered significant using Benjamini-Hochberg FDR (corrected p < 0.05). String analysis [43] was carried out to show the protein interactive network. Figures were plotted in the R environment (v3.1.2) (http://www.r-project.org).