A Comparative Transcriptome Analysis Reveals Physiological Maturation Properties of Mycelia in Pleurotus tuoliensis

Pleurotus tuoliensis is a precious edible fungus with extremely high nutritive and medicinal value. The cultivation period of P. tuoliensis is longer than those of other Pleurotus species, which is mainly due to a longer mycelium physiological maturation period (30–60 days). Currently, the molecular processes underlying physiological maturation of the mycelium remain unclear. We performed a comparative transcriptomic analysis of immature and mature mycelia using RNA-seq. De novo transcriptome assembly resulted in identification of 17,030 unigenes. 451 differentially expressed genes—including those encoding nucleoside diphosphate kinase (NDPK), glycoside hydrolase family proteins, exopolygalacturonase, and versatile peroxidases—were identified. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses revealed that nucleotide synthesis and energy metabolism are highly active during the physiological maturation of mycelia, and genes related to these pathways were significantly upregulated in mature mycelia. NDPK is predicted to be essential for mycelia maturation. Our findings contribute to a comprehensive understanding of mycelia maturation in a commercially important fungal species. Future efforts will focus on the function of NDPK and the mechanism by which it regulates mycelia maturation.


Introduction
Pleurotus tuoliensis, also known as Bailinggu, is one of the most commercially important species of edible mushrooms [1]. P. tuoliensis was initially found in Xinjiang province [2] and is now widely commercially cultivated in China, Japan, and Korea, due to its desirable nutritional and medicinal value [3].
In recent years, P. tuoliensis yield decreased year by year. In 2017, only 70,000 tons were produced, which represents a 77% decrease from 2012 (http://hz.cefa.org.cn/2019/03/18/10498.html). This decrease may be due to this fungus' biological characteristics, such as long cultivation period, a high proportion of deformed fruiting bodies, poor uniformity of fruiting bodies, and strain degeneration [4]. Together, these factors present challenges for mushroom farmers and thus lead to a decline in cultivation enthusiasm. Physiological maturation of P. tuoliensis mycelia is essential for its fruiting. This maturation takes 30-60 days and is responsible for the above-mentioned long cultivation period. In this physiological maturation stage, mycelia accumulate and metabolize nutrients [5]. Cold stimulation of physiologically mature mycelia can result in improved consistency in budding of primordia in culture bottles and ensures a high degree of fruiting body uniformity [6]. Our previous study demonstrated that

Collection of Pleurotus Tuoliensis Mycelia
P. tuoliensis (ACCC 50869) was obtained from the Agricultural Culture Collection of China. Mycelia were grown at 25 • C in plastic bottles containing 820 g of cultivation substrate (40.43% cotton seed hull, 21.56% corncob, 25.16% wheat bran, 9.88% corn flour, 0.99% calcium carbonate, and 1.98% lime, and 65% water content, pH 8. 5-8.8) in the dark for 37 days, after which the substrate was fully covered by mycelia. We randomly selected three cultivation bottles and took the mycelia from cultivation substrate into three corresponding tubes and froze them at −80 • C until RNA extraction (control samples, A). The remaining plastic bottles were maintained at a temperature of 17 • C in the dark to induce physiological maturation of mycelia. After 20 days of induction, CA-asp content in mycelia were determined every two days. After CA-asp content in mycelia reached 0.9 g/L, mycelia reached physiological maturation [7]. We then randomly selected three cultivation bottles and transferred the mycelia to three corresponding tubes. All mycelia samples were stored at −80 • C until RNA extraction (case samples, B). The remaining plastic bottles were transferred to a mushroom house for fruiting management.

Library Preparation and RNA-Seq
Mycelia total RNA was extracted using TRIzol reagent (Life Technologies, New York, NY, USA) according to the manufacturer's protocol. The concentration, quality, and integrity of the extracted RNA samples were determined by agarose gel electrophoresis and a NanoDrop spectrophotometer (Thermo Scientific, USA). Six cDNA libraries were constructed using the Illumina TruSeq Stranded mRNA LT Library Preparation Kit (Illumina, USA) and then paired-end sequenced on an Illumina HiSeq 2000 platform by Shanghai Personal Biotechnology Co. Ltd. (Shanghai, China).

De Novo Transcriptome Assembly and Annotation
Raw reads were assessed for quality using FastQC (https://www.bibsonomy.org/bibtex/ f230a919c34360709aa298734d63dca3) and have been deposited in the National Center for Biotechnology Information (NCBI) database under the the BioProject accession PRJNA551095. Low-quality reads, duplicate sequences and adapter sequences were removed to obtain high-quality clean reads using cutadapt (https://cutadapt.readthedocs.io/en/stable/). Clean paired-end reads were assembled into contigs based on the similarity between overlapping regions using Trinity [13]. The resulting contigs were clustered with TGICL-2.1 [14] with the parameters "-l40-v20". The transcript isoforms with the longest sequence in each Trinity subcomponent were extracted as a representative sequence of the Trinity subcomponent, called Unigene. All unigenes were pooled to yield the final unigene library.
Based on sequence homology, all assembled unigenes were compared to the NCBI non-redundant protein (Nr) and Swiss-Prot databases using the BLASTX program (e-value ≤ 1 × 10 −5 ). The Blast2GO program [15] was used to obtain gene ontology (GO) annotations for the unigenes, and GO-term classification was conducted based on the Nr annotations. Pathway annotations were assigned using the BBH (bi-directional best hit) method of the Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatic Annotation Server (KAAS) online tool [16]. Unigenes were also aligned to the Evolutionary Genealogy of Genes: Non-Supervised Orthologous Groups (eggNOG) database for further functional prediction and classification.

Identification and Functional Analysis of Differentially Expressed Genes (DEGs)
To identify the differentially expressed genes, transcript abundance were first calculated by mapping clean reads to unigenes database using RSEM [17]. DEGs were then identified by R Bio-package DEseq2. [18]. The significance of gene expression differences was assessed using the p-value ≤ 0.05 and |log2 ratio| ≥ 1. These differentially expressed data have been deposited in the Gene Expression Omnibus (GEO) database at the NCBI archives (https://www.ncbi.nlm.nih.gov/geo) under accession GSE135839.
The online analysis tool, the Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) were used to analysis DEGs for GO term and KEGG pathyway enrichment.

Reverse Transcription-Quantitative PCR (RT-qPCR) Analysis
We selected 12 DEGs identified by RNA-seq for validation by RT-qPCR. Total RNA from immature and mature mycelia was extracted using the same procedures as used for RNA-seq. RNA from each sample was reverse transcribed using TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China). Primers for the 12 DEGs were designed in Primer 6.0 software (Premier, Ottawa, Canada) and are shown in Table 1. Reaction mixtures (20 µL) contained 10 µL of 2× SYBR Green Master Mix, 0.4 µL of 10 nM each primer, 2 µL of cDNA template, and 7.2 µL of ddH 2 O. The qPCR reactions were cycled on an ABI Prism 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) with SYBR green fluorescent dye (TaKaRa, Dalian, China). Amplification conditions were: 95 • C for 3 min, followed by 40 cycles of 94 • C for 3 s, 60 • C for 32 s, and 72 • C for 30 s. Each reaction was performed in triplicate using the glyceraldehyde-phosphate dehydrogenase gene (gapdh) as an internal reference. Relative gene expression levels were calculated using the 2 −∆∆Ct method [19].

Growth of P. Tuoliensis Mycelia
The cultivation substrate was fully covered by white fungal mycelia after incubation at 25 • C for 37 days in the dark. These cultivation bottles were then maintained at a temperature of 17 • C in the dark to induce mycelia maturation. The sawdust media surface was gradually covered with plump and elastic hyphae (Figure 1). At the 35th day, CA-asp content in mycelia was determined to be 0.982 g/L, which means the physiological maturation of mycelia. Primordia were initiated within 7 days under low temperature and light stimulation conditions and within 15 days, developed reproductive fruiting bodies that could be harvested. To investigate the molecular basis of physiological maturation of mycelia, we collected immature mycelia and mature mycelia for RNA-seq analysis.

Growth of P. Tuoliensis Mycelia
The cultivation substrate was fully covered by white fungal mycelia after incubation at 25 °C for 37 days in the dark. These cultivation bottles were then maintained at a temperature of 17 °C in the dark to induce mycelia maturation. The sawdust media surface was gradually covered with plump and elastic hyphae (Figure 1). At the 35th day, CA-asp content in mycelia was determined to be 0.982 g/L, which means the physiological maturation of mycelia. Primordia were initiated within 7 days under low temperature and light stimulation conditions and within 15 days, developed reproductive fruiting bodies that could be harvested. To investigate the molecular basis of physiological maturation of mycelia, we collected immature mycelia and mature mycelia for RNA-seq analysis.

Sequencing and Gene Annotation
We performed paired-end sequencing of six cDNA libraries (three biological replicates each for immature and mature mycelia). We obtained an average of 45,688,081 raw reads per library. After removing adapter sequences and low-quality reads, clean reads represented approximately 99.53% of the raw reads. Clean reads were combined into a single dataset for de novo assembly. A total of 60,330 transcripts, including sequence isoforms, were obtained. Transcripts displayed a mean length of 1897 bp. To reduce redundancy, only the longest sequence isoform within each subcomponent

Sequencing and Gene Annotation
We performed paired-end sequencing of six cDNA libraries (three biological replicates each for immature and mature mycelia). We obtained an average of 45,688,081 raw reads per library. After removing adapter sequences and low-quality reads, clean reads represented approximately 99.53% of the raw reads. Clean reads were combined into a single dataset for de novo assembly. A total of 60,330 transcripts, including sequence isoforms, were obtained. Transcripts displayed a mean length of 1897 bp. To reduce redundancy, only the longest sequence isoform within each subcomponent was selected. In total, we obtained 17,030 unigenes representing unique transcripts for further analysis ( Figure 2).
Genes 2019, 10, x FOR PEER REVIEW 5 of 15 was selected. In total, we obtained 17,030 unigenes representing unique transcripts for further analysis ( Figure 2). This unigene dataset displayed a total size of 21,120,061 bp, an average length of 1240 bp, and an N50 of 1891 and served as the mRNA transcriptome for mycelia ( Table 2).  This unigene dataset displayed a total size of 21,120,061 bp, an average length of 1240 bp, and an N50 of 1891 and served as the mRNA transcriptome for mycelia (Table 2).
In total, 14,300 and 7115 unigenes were mapped using the eggNOG (Table S1) and Swiss-Prot (Table S2) databases, respectively. In addition, 13,611 (79.92%) unigenes matched known proteins in NR (NCBI non-redundant protein sequences) database (Table S3), and 95% of these were homologous with proteins from P. ostreatus, a species evolutionarily related to P. tuoliensis (Table S4). Additionally, 8,804 (51.70%) unigenes assigned to more than one GO term were classified into 67 functional groups from the Gene Ontology database (Table S5). Using KAAS (KEGG Automatic Annotation Server) with P. ostreatus as the reference, a total of 2,942 unigenes were mapped to 35 level-2 pathways (Table S6). Furthermore, a total of 2603 unigenes were annotated in all of these databases. a Percentage of the number of clean reads in the number of total raw reads. b All sequences are arranged from long to short sequence, and the sequence lengths are sequentially added in this order, When the added length reaches 50% of the total length of the sequence, the length of the last sequence is recorded as N50 (bp).

Identification and Functional Enrichment Analysis of Differentially Expressed Genes (DEGs)
To identify significantly differentially expressed genes (DEGs), we compared the normalized read counts (FPKM values) of the immature and mature mycelia samples. This comparative analysis revealed 451 significantly DEGs (p-value ≤ 0.05, |log 2 ratio| ≥ 1), which included 206 upregulated and 245 downregulated genes ( Figure 3). The analyzed data of differential expression analysis for each unigene are shown in Table S7.
We next performed GO analysis on the DEGs to gain insight into their biological functions. A total of 197 DEGs were enriched. Enrichment result of all DEGs is shown in Table S8. Figure 4 shows the top 32 enriched GO terms assigned to our DEGs. The cellular metabolic process (GO:0044237) annotation was enriched with most DEGs (False discovery rate (FDR) > 0.05). Oxidation-reduction process (GO:0055114) and oxidoreductase activity (GO:0016491) were also enriched (FDR > 0.05), indicating that oxidation-reduction processes may be highly active in mature mycelia. Other enriched GO terms were mainly involved in the synthesis and processing of genetic information, including rRNA modification (GO:0000154) (FDR > 0.05), mRNA pseudouridine synthesis (GO:1990481) (FDR < 0.05), snRNA metabolic process (GO:0016073) (FDR > 0.05), and deoxyribonucleoside diphosphate metabolic process (GO:0009186) (FDR > 0.05). Within the molecular function category, terms associated with carbohydrate metabolism, such as hydrolase activity acting on O-glycosyl compounds (GO:0004553) (FDR > 0.05) and hydrolase activity acting on glycosyl bonds (GO:0016798) (FDR > 0.05), were enriched. In general, GO terms associated with genetic information synthesis and processing were largely enriched.
To identify significantly differentially expressed genes (DEGs), we compared the normalized read counts (FPKM values) of the immature and mature mycelia samples. This comparative analysis revealed 451 significantly DEGs (p-value ≤ 0.05, |log2 ratio| ≥ 1), which included 206 upregulated and 245 downregulated genes ( Figure 3). The analyzed data of differential expression analysis for each unigene are shown in Table S7.  Table S8. Figure 4 shows the top 32 enriched GO terms assigned to our DEGs. The cellular metabolic process (GO:0044237) annotation was enriched with most DEGs (False discovery rate (FDR) > 0.05). Oxidation-reduction process (GO:0055114) and oxidoreductase activity (GO:0016491) were also enriched (FDR > 0.05),  To better understand the functions and interactions of our DEGs, all DEGs underwent a pathway-based analysis by mapping against the KEGG database. Pathway enrichment analysis demonstrated that only 35 pathways were enriched, which included 45 DEGs (Table S9). We found that these pathways were mainly associated with nucleotide metabolism, DNA replication and repair, and energy metabolism, such as pyrimidine metabolism (ko00240), nucleotide excision repair (ko03420), pyruvate metabolism (ko00620), and the pentose phosphate pathway (ko00030). The To better understand the functions and interactions of our DEGs, all DEGs underwent a pathway-based analysis by mapping against the KEGG database. Pathway enrichment analysis demonstrated that only 35 pathways were enriched, which included 45 DEGs (Table S9). We found that these pathways were mainly associated with nucleotide metabolism, DNA replication and repair, and energy metabolism, such as pyrimidine metabolism (ko00240), nucleotide excision repair (ko03420), pyruvate metabolism (ko00620), and the pentose phosphate pathway (ko00030). The results of our KEGG analysis corresponded with those obtained by our GO enrichment analysis. The 20 enriched pathways with FDR < 1 are shown in Figure 5.
In addition, the total annotation results of Unigenes and DEGs are listed in Table S10.

Identification of DEGs Related to Physiological Maturation of Mycelia
To identify genes related to physiological maturation of mycelia, we analyzed significantly upregulated DEGs in mature mycelia. Overall, these upregulated DEGs are mainly involved in the synthesis of genetic materials and energy metabolism, and include unigenes encoding glutamate decarboxylase Gad1 (DN579_c0_g1), oxalate decarboxylase OxdC (DN8634_c0_g1), pseudouridylate synthase Pus (DN2821_c0_g2), DNA repair and recombination protein RHM52 (DN3170_c0_g1), and nucleoside diphosphate kinase Ndpk (DN4800_c0_g2), demonstrating that this complex developmental program requires not only fast mycelial vegetative growth, but also more energy to accumulate nutrients and thus nurture developing fruit bodies during the reproductive stage. Versatile peroxidases (Vps) are key lignin-degrading enzymes, as they possess both manganese and lignin peroxidase activities. We observed upregulation of Vp-encoding genes (DN9023_c0_g1, DN9023_c0_g3, and DN8870_c0_g2) in mature mycelia, consistent with nutrient accumulation. Another important class of DEGs identified in our mature mycelia were those encoding post-modification enzymes, such as cytochrome P450 monooxygenase and methyltransferases, which are important for the synthesis of secondary metabolites. The expression of genes encoding DEG number indicates the number of differential genes annotated in the pathway, indicated by a dot, the larger the dot, the more differential genes are enriched. The color gradient indicates the corrected p-value by Benjamini-Hochberg FDR method. Only pathways with FDR < 1 were shown.

Identification of DEGs Related to Physiological Maturation of Mycelia
To identify genes related to physiological maturation of mycelia, we analyzed significantly upregulated DEGs in mature mycelia. Overall, these upregulated DEGs are mainly involved in the synthesis of genetic materials and energy metabolism, and include unigenes encoding glutamate decarboxylase Gad1 (DN579_c0_g1), oxalate decarboxylase OxdC (DN8634_c0_g1), pseudouridylate synthase Pus (DN2821_c0_g2), DNA repair and recombination protein RHM52 (DN3170_c0_g1), and nucleoside diphosphate kinase Ndpk (DN4800_c0_g2), demonstrating that this complex developmental program requires not only fast mycelial vegetative growth, but also more energy to accumulate nutrients and thus nurture developing fruit bodies during the reproductive stage. Versatile peroxidases (Vps) are key lignin-degrading enzymes, as they possess both manganese and lignin peroxidase activities. We observed upregulation of Vp-encoding genes (DN9023_c0_g1, DN9023_c0_g3, and DN8870_c0_g2) in mature mycelia, consistent with nutrient accumulation. Another important class of DEGs identified in our mature mycelia were those encoding post-modification enzymes, such as cytochrome P450 monooxygenase and methyltransferases, which are important for the synthesis of secondary metabolites. The expression of genes encoding cytochrome P450s monooxygenase Psih (DN1132_c0_g1, DN1330_c0_g1, DN5495_c0_g1, DN7493_c0_g1, DN8894_c1_g1, DN8894_c1_g2, and DN9463_c0_g1) and a methyltransferase (DN4644_c0_g2) were found to be significantly upregulated in mature mycelia. Finally, upregulated fruiting body protein SC3 (DN5369_c0_g1) was enriched in mature mycelia, probably indicating the initiation and formation of fruiting bodies (Figure 6).

Discussion
The physiological maturation stage of mycelia, which lasts approximately 30-60 days, is critical for the growth and development of the P. tuoliensis fruiting body. Only mature mycelia are subjected to low temperature and light stimulation, primordia and fruiting body can be finely developed. Transcriptomic studies on mature mycelia to cold stimulation and fruiting body formation for P. tuoliensis have been performed by Fu et al. (2016Fu et al. ( , 2017 to reveal physiological metabolic properties underlying primordia formation and fruiting body development [6,12]. Transcriptomic analysis on cold stimulation of P. tuoliensis mycelia indicated that functional groups of differentially expressed unigenes associated with cell wall and membrane stabilization, soluble sugars, and protein biosynthesis and metabolism pathways, and calcium signaling and mitogen-activated protein kinases (MAPK) pathways play a vital role in Bailinggu's response to cold stimulation. While transcriptomic analysis on fruiting body development identified the stage-specific and differentially expressed unigenes (DEGs) involved in morphogenesis, primary carbohydrate metabolism, cold stimulation, and blue-light response, predicting that these unigenes might help P. tuoliensis adapt to genetic and environmental factors that influence fructification. However, no transcriptome analyses have been performed to study the molecular mechanisms underlying physiological maturation of mycelia. In this study, we utilized a comparative transcriptomic approach to explore differential gene expression and molecular processes during the physiological maturation of mycelia.
Our comparative transcriptome analysis identified 451 differentially expressed genes (DEGs) between immature and mature mycelia. Although the number of DEGs is small considering ~17000 Strain of Pleurotus tuoliensis

Discussion
The physiological maturation stage of mycelia, which lasts approximately 30-60 days, is critical for the growth and development of the P. tuoliensis fruiting body. Only mature mycelia are subjected to low temperature and light stimulation, primordia and fruiting body can be finely developed. Transcriptomic studies on mature mycelia to cold stimulation and fruiting body formation for P. tuoliensis have been performed by Fu et al. (2016Fu et al. ( , 2017 to reveal physiological metabolic properties underlying primordia formation and fruiting body development [6,12]. Transcriptomic analysis on cold stimulation of P. tuoliensis mycelia indicated that functional groups of differentially expressed unigenes associated with cell wall and membrane stabilization, soluble sugars, and protein biosynthesis and metabolism pathways, and calcium signaling and mitogen-activated protein kinases (MAPK) pathways play a vital role in Bailinggu's response to cold stimulation. While transcriptomic analysis on fruiting body development identified the stage-specific and differentially expressed unigenes (DEGs) involved in morphogenesis, primary carbohydrate metabolism, cold stimulation, and blue-light response, predicting that these unigenes might help P. tuoliensis adapt to genetic and environmental factors that influence fructification. However, no transcriptome analyses have been performed to study the molecular mechanisms underlying physiological maturation of mycelia. In this study, we utilized a comparative transcriptomic approach to explore differential gene expression and molecular processes during the physiological maturation of mycelia.
Our comparative transcriptome analysis identified 451 differentially expressed genes (DEGs) between immature and mature mycelia. Although the number of DEGs is small considering~17000 unigenes, this result is reasonable. Even if P. tuoliensis mycelia have undergone a physiological maturation process, their morphology does not display a huge change relative to the process of primordia formation or fruiting body development. As shown in GO and KEGG enrichment results of DEGs, DEGs were mainly enriched for nucleotide synthesis and energy metabolism functions. Which are consistent with the fast growth and nutrient accumulation that are characteristic of mycelium maturation [21]. Metabolome analysis of physiological maturation of P. tuoliensis mycelia also revealed similar metabolism phenomenon. Metabolites enriched on pyrimidine synthesis were significantly abundant in mature P. tuoliensis mycelia, such as CMP and cytosine, which are indicative of the high demand for rapid nucleotide synthesis during physiological maturation [7].
Among our DEGs, the upregulation expression of NDPK is interesting. Many studies have demonstrated that NDPK regulates various biological processes and signal transduction pathways in addition to playing a fundamental role in maintaining the cellular balance of NDP and NTP [22,23]. In Neurospora crassa, Ndpk is associated with hyphal development and photomorphogenesis [24]. Tang et al. (2016) reported that Ndpk was significantly upregulated in mature Lentinus edodes mycelia, concluding that this enzyme may plays an important role in light-induced mycelial brown-film formation, an indicator of physiological mycelial maturation in L. edodes [21]. Thus, upregulation of NDPK in mature P. tuoliensis mycelia is likely associated with physiological mycelial maturation. The specific function of Ndpk in this process is worthy of further study.
Growth and development of P. tuoliensis relies on nutrients produced by decomposition of the cultivation substrate. CAZymes are involved in the hydrolysis of plant cell wall polysaccharides and play an important role in substrate degradation processes [25]. We identified a total of 21 CAZymes, most of which belong to glycoside hydrolase (GH) families, which play important roles in lignocellulose decomposition [26]. The GH families GH13 (2 genes), GH15 (1 gene), GH28 (3 genes), GH5 (2 genes), and GH1 (1 gene) were upregulated in our mature mycelia. GH5 is one of the largest GH families, whose members act on a wide range of substrates, including β-1,3-glucans, β-1,4-glucans in cellulose, and β-1,4-mannans in hemicelluloses [27]. GH13 is a family of starch-debranching enzymes that hydrolyze the α-1,6-glucosidic bond in starch, and GH28 members play critical roles in pectin degradation [28]. The function of GH15 is not known. Lignin degradation involves three types of degrading enzymes: manganese peroxidases (MnPs) [29], versatile peroxidases (VPs) [30], and laccases [31]. We identified three upregulated VP-encoding genes in mature mycelia. These upregulated genes are probably related with lignin degradation.
Functional enrichment analyses of our DEGs also demonstrated that oxidation-reduction activities and pathways are very active during the physiological maturation of mycelia. Cytochrome P450 monooxygenases psiH, which catalyze the biooxidation of various substrates through activation of molecular oxygen and play important roles in metabolic processes and stress responses [32], displayed increased expression in mature mycelia. Similar expression patterns for these genes were observed by qPCR analysis. Identification of upregulated PSIH demonstrated that mycelium maturation requires increased oxygen and produces excess water, which often leads to submergence and hypoxic stress.
We also observed upregulated fruiting body protein SC3 in mature mycelia, which was reported to be involved in fruiting body development of edible mushroom [33,34]. The high abundance of this protein probably indicates the initiation and formation of fruiting bodies. Fu et al. (2017) found that unigenes encoding fruiting body hydrophobin 1 were significantly upregulated in primordia and fruiting body for P. tuoliensis. The unigene encoding hydrophobin 2 was upregulated in primordia while downregulated in the fruiting bodies. However, they did not find unigenes that specifically expressed during the initiation of fruiting [12].
In summary, our comparative transcriptomic analysis of immature and mature mycelia revealed gene expression differences indicative of changing metabolism during physiological mycelium maturation and identified a potential functional group of genes responsible for regulating the physiological maturation in P. tuoliensis mycelia. GO and KEGG enrichment analysis demonstrated that metabolic pathways related to genetic information and energy metabolism are extremely active during the physiological maturation of mycelia. Genes related to critical oxidation-reduction processes are upregulated in mature mycelia. CAZymes involved in the hydrolysis of plant cell wall polysaccharides were differentially expressed in mature mycelia compared with immature mycelia and the differential expression of these genes may relate with substrate degradation utilization. In addition, we predict that NDPK may be an important enzyme for participating in mycelia physiological maturation. Our data enables an improved understanding of physiological maturation in P. tuoliensis mycelia, which lay an important foundation for revealing mycelia physiological maturation properties in edible mushrooms.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/9/703/s1, Table S1: eggNOG annotation statistics of unigenes; Table S2: unigenes annotation information in Swiss-Prot databases; Table S3: unigenes matched known proteins in NR database; Table S4: Nr annotation statistics of unigenes; Table S5: GO enrichment terms of unigenes; Table S6: KEGG pathways enriched with unigenes; Table S7: the total data of differential expression analysis for each unigene; Table S8: GO enrichment terms of differential expressed genes; Table S9: KEGG enrichment pathways of differential expressed genes; Table S10: the total annotation of unigenes and differential expressed genes.