Comparative transcriptome analysis of dikaryotic mycelia and mature fruiting bodies in the edible mushroom Lentinula edodes

Lentinula edodes is a popular cultivated edible mushroom with high nutritional and medicinal value. To understand the regulation of gene expression in the dikaryotic mycelium and mature fruiting body in the commercially important Korean L. edodes strain, we first performed comparative transcriptomic analysis, using Illumina HiSeq platform. De novo assembly of these sequences revealed 11,675 representative transcripts in two different stages of L. edodes. A total of 9,092 unigenes were annotated and subjected to Gene Ontology, EuKaryotic Orthologous Groups, and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Gene expression analysis revealed that 2,080 genes were differentially expressed, with 1,503 and 577 upregulated in the mycelium and a mature fruiting body, respectively. Analysis of 18 KEGG categories indicated that fruiting body-specific transcripts were significantly enriched in ‘replication and repair’ and ‘transcription’ pathways, which are important for premeiotic replication, karyogamy, and meiosis during maturation. We also searched for fruiting body-specific proteins such as aspartic protease, gamma-glutamyl transpeptidase, and cyclohexanone monooxygenase, which are involved in fruiting body maturation and isolation of functional substances. These transcriptomes will be useful in elucidating the molecular mechanisms of mature fruiting body development and beneficial properties, and contribute to the characterization of novel genes in L. edodes.

different stages, e.g., mycelia, brown mycelial film, primordium, and the fruiting body 18,19 , to determine the stage at which optimum activity or the maximum yield of the desired substance occurs. Therefore, determining genetic changes compared with the parental strains at different developmental stages using high-throughput sequencing technology (next-generation sequencing technology) 20 will be necessary to generate improved strains and varieties of mushrooms by genetic modification.
To analyze differentially expressed genes (DEGs) between two major developmental stages (i.e., mycelium and fruiting body), high-throughput sequencing has been used for genome and transcriptome analyses of basidiomycota mushrooms, including Agrocybe aegerita, Auricularia polytricha, Cordyceps militaris, and Ganoderma lucidum. Fruiting body formation is regulated at the transcriptional level and results in morphological changes during development and stage-specific expression to produce beneficial properties [21][22][23][24] . To date, there have been no transcription-level, genome-wide comparative studies of the development stages of mycelia and mature fruiting bodies of L. edodes commercial strains and their relation to yield. Two monokaryotic strains of L. edodes have been sequenced and annotated by high-throughput sequencing 25,26 and extensive expressed gene sequence data of a Chinese strain was obtained by deep Solexa sequencing 27 . Differences in expression between the dikaryotic mycelium and primordium have also been investigated by serial analysis of gene expression 28 . Comparative transcriptome analysis has been performed to assess fruiting bodies at different stages of development 19,29 . In addition, transcriptome analyses have been conducted in light-induced brown film, different growth stages of the fruiting body, and a post-harvested fruiting body 16,17,19 , and post genome-wide studies have been performed to understand the complex molecular mechanisms underlying gene regulation in L. edodes. However, despite the great importance of investigating the specific stages for the optimal production of beneficial nutritional and medicinal properties and the generation of improved L. edodes cultivars, the transcriptomes of dikaryotic mycelia and mature fruiting bodies of this species remain poorly understood. In particular, there have been no reports of comparative transcriptomic analysis of a Korean L. edodes cultivar using the NGS technique, despite active culturing and consumption of great quantities of edible L. edodes mushrooms in other East Asian countries. In this study, we performed the first comparative transcriptome analysis of two developmental stages, the mycelium and mature fruiting body.
In this study, we performed the first comparative transcriptome analysis of two developmental stages, the mycelium and a mature fruiting body, of L. edodes using Illumina sequencing technology to explore the transcriptome during maturation of the fruiting body at the genome level. Approximately 30 million reads were obtained from each sample and mapped to the L. edodes genome. A total of 11,675 unigenes were identified by de novo assembly, and analyses of the DEGs revealed genes involved in fruiting body maturation of L. edodes. A functionally annotated and classified transcriptome dataset will provide valuable genomic information for further studies in L. edodes. The identification of differentially expressed genes (DEGs) involved in mycelium and fruiting body development will improve our understanding of the maturation and stage-specific expression of functional properties of this edible mushroom, and provide valuable genetic information for further studies of the molecular mechanism aiming to improve productivity.

Results and Discussion
Transcriptome sequencing of L. edodes and de novo assembly. To identify transcriptomic changes in two developmental stages of the commercially important Korean L. edodes strain Sanjo701 (SJ701) 3,4 , high-throughput sequencing in RNA samples extracted from L. edodes dikaryotic mycelia and mature fruiting bodies was the first performed using Illumina paired-end sequencing technology according to the transcriptome analysis workflow (Fig. 1).
In total, 27,646,278 and 30,109,768 clean paired-end sequence reads were obtained from the mycelium and mature fruiting body, respectively, with a Q20 percentage >72.2%, using the DynamicTrim and LengthSort programs 30 . Trimming data resulted in an average read length of 82.8 bp across all samples, with a minimum read length of 25 bp, using de novo transcriptome assemblers. A total of 57,756,046 high-quality reads (4,780,938,320 bp) were obtained after combining the mycelium and mature fruiting body libraries of the dikaryon L. edodes strain SJ701 (ASI no. 3305) ( Table 1). These reads were then assembled into 32,001 extended transcripts (k-mer = 51) with a mean length of 2,152 bp and N50 length of 2,900 bp ( Table 2). The transcript length ranged from 200 to >4,500 bp ( Supplementary Fig. S1). Finally, 11,675 loci were predicted from the extended transcripts, and the average length of the assembled loci was 1,702 bp with an N50 length of 2,443 bp ( Table 2).
The genome size of L. edodes ranged from 41.8 to 46.1 Mb 25,26 . As a result, approximately 4.78 Gbp of transcriptome sequence data represented at least 103-fold coverage. A reference-based assembly should reconstitute full-length transcripts with more than 10-fold coverage, and higher coverage is considered more challenging 31 . Therefore, the high coverage in this study will allow for more accurate assembly of the transcriptome. According to a recent report, 13,028 genes in L. edodes have been predicted by whole genome de novo sequencing and genome annotation 26 . Despite identifying approximately 89.6% (11,675) of these previously predicted genes (13,028), our results will allow for further annotation and characterization of the dikaryotic L. edodes strain.
Functional annotation and classification of L. edodes transcripts. Genome annotation and functional classification of 9,092 unigenes (77.9% of 11,675) were determined from the six datasets for the L. edodes SJ701 strain. All unigenes were aligned against sequences from each gene database using the BLASTx algorithm (e-value ≤ 1e −10 ). A total of 8,980 (98.8%) unigenes were found to match known proteins in the NR database, as well as the UniProtKB fungi (8,190,90   To assess the enrichment of DEGs between L. edodes mycelia and mature fruiting bodies, we performed Gene Ontology (GO), Eukaryotic Orthologous Groups (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. The transcripts were functionally classified as follows. First, GO enrichment analysis is a bioinformatics tool used to annotate genes and gene products to terms in specific ontologies. GO terms are composed of three primary ontologies: biological process, cellular component, and molecular function 33 . To classify the functions of the annotated unigenes based on those in other fungi, 5,834 of the 9,092 unigenes were classified into the three primary GO categories and 18 functional subcategories, with 6, 6, and 6 annotated to biological process, cellular component, and molecular function, respectively, at an ontology depth score of 1 34 (Supplementary  Table S3). To investigate further, the level of differentiated depth was used for GO classification, unlike the depth of other reported results in L. edodes 19,27 . A total of 5,834 unigenes were assigned to at least one GO term. Cellular process (2,698 unigenes; 45.8%) and metabolic process (2,415; 41.0%) were the most common terms in the various biological process subcategories, cell (2,479; 75.3%) and membrane (749; 22.8%) were the most common terms in the cellular component category, and catalytic activity was the most common term in the molecular function category (Supplementary Table S3). These five GO terms (88.8% of all GO terms) represented the majority of the subcategories in the GO classification. Furthermore, the major processes represented by the L. edodes GO terms were consistent with previously reported results of comparisons between mycelia and fruiting bodies in A. aegerita, A. polytricha, Hypsizygus marmoreus, and Pleurotus tuoliensis 21,22,35,36 . In contrast, GO functional classification in two different studies 19,27 , which simply characterized the mycelia of a single strain and compared gene expression based on three different stages of fruiting body growth in L. edodes, showed significantly different results, despite the application of similar transcriptomic analyses. Thus, the current study is the first to demonstrate that L. edodes gene expression is regulated very differently depending on the developmental stage.
Second, to further analyze the annotations and transcriptome, the KOG database 36,37 was used to assign the unigenes to 25 classifications for the first time in L. edodes. A total of 4,904 unigenes were assigned to at least one term in the 25 KOG functional categories (Supplementary Table S4). The top five KOG categories included 'posttranslational modification, protein turnover, chaperones' (O; 429, 8.7%), 'signal transduction mechanisms' (T; 358, 7.3%), 'secondary metabolites biosynthesis, transport, and catabolism' (Q; 299, 6.1%), 'translation, ribosomal structure, and biogenesis' (J; 292, 6.0%), 'nuclear structure' (X; 292, 6.0%) with the exception of the largest KOG group, and 'general function prediction only' (R; 788, 16.1%) (Supplementary Table S4). The top 10 KOG categories were matched to at least seven of the clusters of orthologous groups previously reported in L. edodes with respect to the Basidiomycota A. aegerita, A. polytricha, H. marmoreus, and P. tuoliensis 19,21,22,27,35,36 . However, the ranking of annotated categories differed. Interestingly, KOG classifications showed a higher similarity among other studies of L. edodes if several comparable unigenes were used for the clusters of orthologous group or KOG annotations.
Finally, to identify biological pathways in L. edodes, the 11,675 transcript sequences were searched in the KEGG database, and 1,717 unigenes were assigned to KEGG pathways, which included 4 main categories with 18 subcategories (98 pathways). The pathways consisted of metabolism-(875 unigenes, 51.0%), genetic information processing-(586, 34.1%), cellular processes-(166, 9.7%), and environmental information processing-related (85, 5.0%) functional categories, which are associated with the maintenance of fundamental biological pathways (Supplementary Table S5). Among the major categories, metabolic pathways, especially carbohydrate and amino acid metabolism, were active in L. edodes, similar to other mushrooms including A. aegerita, A. polytricha, and G. lucidum, in which comparative transcriptome analyses between mycelia and mature fruiting bodies have been performed 21,22,24 . Surprisingly, active subcategories of the KEGG annotated results based on transcriptome analysis data from the mycelium and mature fruiting body were more similar to those of the same stage in other mushrooms 21,22,24 than to different stages of L. edodes 19,27 . The KEGG analysis results were similar to those of the GO functional analysis in this study, suggesting that KEGG annotation and analysis are valuable for exploring DEGs during different developmental stages. This observation has not been previously described in transcriptome analyses of the mycelial stage and the three stages of fruiting body growth in L. edodes. Therefore, the GO, KOG, and KEGG analyses in this study will be useful for further expression profiling and determination of DEGs, and these analyses are significant compared with previously reported results in L. edodes.

Analysis of DEGs involved in the different developmental stages of L. edodes. To identify DEGs
between the dikaryotic mycelium and mature fruiting body of L. edodes, the expression level of each unigene from the two developmental stages was firstly calculated using the DESeq package 38 and Another Multidimensional Analysis Package 39 , setting a threshold value of log 2 ratio ≥ 1 and false discovery rate ≤ 0.01. The MA (scatter plot of the log 2 fold change vs. the mean count) plot of pairwise comparisons is shown in Fig. 3A, in which red and green dots indicate the detected DEGs. Of 11,675 unigenes, 2,080 were significantly differentially expressed between the two stages in L. edodes, with 1,503 and 577 unigenes upregulated in the mycelium and mature fruiting body, respectively (Fig. 3B), and 1,211 unigenes expressed in both stages (Fig. 4). In total, 532 and 337 unigenes were specifically expressed in the mycelium and mature fruiting body, respectively. Overall, 9,595 unigenes showed significant changes in expression between the dikaryotic mycelium and mature fruiting body (Fig. 4). This expression profile is inconsistent with previous reports of DEGs detected in the dikaryotic mycelium and primordium by serial analysis of gene expression (SAGE) 28 , fruiting body formation phases by RNA sequencing (RNA-Seq) 19 , and mycelia and/or fruiting bodies cultured under different conditions 16,25 due to large differences in gene regulation. The expression conditions were also inconsistent. Similar numbers of stage-specific DEGs of the mycelia and primordia were found in a gene profiling analysis using SAGE 28 . Different numbers of DEGs in three different stages during fruiting body growth were observed in a pairwise comparison of GO annotations using high-throughput RNA-Seq; in particular, genes were more highly upregulated in the mature fruiting body than in the primordium 19 . However, no genetic information related to DEGs of the mycelium and fruiting body has yet been reported. Therefore, our results necessitate further exploration of L. edodes DEGs in future studies. To investigate the function of the DEGs, the identified genes were annotated by GO, KOG, and KEGG. GO enrichment analysis provided functional GO terms that were substantially enriched in DEGs compared with the genome background. The results indicate that the DEGs were connected to interesting biological functions. To obtain detailed genetic information on the DEGs using the GO database, we classified the three primary GO categories into 86 functional subcategories, and 40, 17, and 29 unigenes were annotated in the biological process, cellular component, and molecular function categories, respectively, at an ontology depth score >3 (Fig. 5). Interestingly, the number of unigenes with GO annotations upregulated in the maturing fruiting body (407) was significantly less than the number of downregulated unigenes (1597) as transcriptomic studies in the same stage of C. militaris 23 . The number of DEGs annotated to the three primary GO terms was greater in the mycelium than in the mature fruiting body. In particular, the number of DEGs involved in biological processes was increased more than six-fold in mycelia (735) compared with the mature fruiting body (116) (Supplementary Table S6). These results suggest that the gene expression profile varies widely during mycelial growth compared with fruiting body maturation in L. edodes. Because no reports of GO enrichment analyses to compare the mycelium and fruiting body during any growth stage have been published, our results are novel and will be useful in the future exploration of DEGs in L. edodes.
GO enrichment analysis revealed that the cellular component 'intracellular membrane-bounded organelle' (248; 12.4% of 2,004 DEGs) and 'cytoplasm subcategories' (239; 11.9%) were most strongly enriched during both developmental stages. Although the terms 'protein metabolic process' (120; 7.5% of mycelial DEGs), 'small molecule metabolic process' (111; 7.0%), and 'cellular aromatic compound metabolic process' (90; 5.6%) in the biological process category were over-represented among the mycelial DEGs, no DEGs in mature fruiting bodies belonged to these three subcategories. 'Peptidase activity' (35 DEGs), 'transferase activity, transferring phosphorus-containing groups' (28), and 'coenzyme binding' (20) in the molecular process GO category were represented among the mycelial DEGs; however, no DEGs in the mature fruiting body were annotated to these subcategories. In contrast, DEGs involved in 'hydrolase activity, acting on acid anhydrides' (18), 'membrane-bounded vesicle' (5), 'mitochondrial membrane' (4), 'response to oxidative stress' (3), and 'response to starvation' (3) were expressed only in the mature fruiting body, according to the GO enrichment analysis (Fig. 5 and Supplementary Table S6). Hydrolase activity is important in cellular and subcellular movement, such as catalyzing transmembrane movement of substances during development of the fungal fruiting body 40,41 . Mitochondria also play a critical role in cell differentiation 42 , and mitochondrial processing enzymes are regulated during formation of the fruiting body 43 . Furthermore, a higher level of oxidative stress and/or nutritional starvation trigger mushroom development in various basidiomycota 35,[44][45][46] . The results of this study suggest that mitochondrial processing, cell remodeling, oxidative stress, and starvation are more active in fruiting body development of the basidiomycete L. edodes.
Enrichment analysis of KOG functional categories was performed to further investigate the difference in classification depending on the developmental stage ( Fig. 6 and Supplementary Table S7). This result has not been described previously. The number of upregulated DEGs was smaller in the maturing fruiting body (305) than in the mycelium (925), similar to the results of the GO enrichment analysis. The top five KOG categories, 'general   for the categories 'replication, recombination, and repair' and 'chromatin structure' . Development of the fruiting body is preceded by premeiotic replication, karyogamy, and meiosis, and thus genes related to these processes are overexpressed during the maturation stages 18,47,48 . These results suggest that genes involved in DNA replication, recombination, repair, chromatin structure, and dynamics are overexpressed during the fruiting body stage.
To further evaluate the functional pathways involved in the fruiting body, all DEGs were matched to the KEGG database ( Fig. 7 and Supplementary Table S8). KEGG sub-classification analysis showed that a greater number of DEGs were involved in carbohydrate and amino acid metabolism and translation, which represent pathways involved in the maintenance of basic biological processes during both stages 22 . Among the KEGG sub-classifications, the number of upregulated DEGs was decreased in the maturing fruiting body compared with the mycelium, consistent with the GO and KOG enrichment analyses. The majority of the sub-classifications, including 'nucleotide metabolism' , 'metabolism of other amino acids' , 'glycan biosynthesis and metabolism' , 'metabolism of cofactors and vitamins' , 'energy metabolism' , 'lipid metabolism' , 'folding, sorting, and degradation' , 'transport and catabolism' , and 'signal transduction' , among all 18 KEGG sub-classifications including the above 3 sub-classifications were over-represented in the mycelium stage (Fig. 7). On the other hand, there was little difference in the number of DEGs in the 'metabolism of terpenoids and polyketides' , 'biosynthesis of other secondary metabolites' , 'membrane transport' , and 'cell growth and death' sub-classifications. Moreover, the 'replication and repair' and 'transcription' sub-classifications were uniquely enriched in the mature fruiting body stage (Fig. 7).
In addition, KEGG pathway analysis was used to compare the differences between the DEGs of both stages to identify a more definitive biological pathway represented among the sub-classifications. Nearly all pathways identified in the KEGG enrichment analysis were more associated with the DEGs in the mycelium than those in the mature fruiting body. Only 6 of the 86 pathways showed a higher number of DEGs in the mature fruiting body compared with the mycelium. In particular, genes involved in 'replication and repair' , as well as the 'mismatch repair' , 'homologous recombination' , 'non-homologous end-joining' , and 'spliceosome of transcription' , were more active in the mature fruiting body compared with the mycelium, similar to the results of the KOG enrichment analysis (Supplementary Table S8). These results are consistent with those described in previous reports, in which genes participating in DNA replication, recombination, repair, chromatin structure, and the associated dynamics were found to be important in development of the fruiting body. As reported previously, spliceosome pathways were was also enriched in the fruiting body of C. militaris 23 . However, the findings of the present study are distinct from previous transcriptome analyses conducted in other mushrooms comparing the mycelium and fruiting body [21][22][23][24] . Gene expression patterns vary widely during mycelial growth compared with fruiting body maturation in L. edodes. Hence, it is possible that downregulation of gene expression is required for mature fruiting bodies in the L. edodes genome, and information regarding fruiting body-related genes is currently limited in the available databases. Special genes involved in mature fruiting body development in L. edodes. Analysis of the DEGs provides valuable information regarding maturation of the fruiting body, which is initiated by the presence of a stimulus that activates specific genes that enable development of the mycelium into a fruiting body 18,[21][22][23][24] . Analysis of the DEGs allowed us to search for marked differences in the gene expression profiles of the mycelium and fruiting body in L. edodes. The log 2 fold change values between the two stages ranged from −11.3 to 11.7, and DEGs showing a greater than twofold change between the two stages are presented in Supplementary Table S8. Overall, 577 unigenes were upregulated in the fruiting body, including 240 co-expressed and 337 fruiting body-specific unigenes (log 2 ratio ≥ 1) (Fig. 4). Of the 337 stage-specific unigenes, 53 unigenes were exclusively overexpressed in the mature fruiting body of L. edodes. The overexpressed fruiting body-specific genes were selected based on read counts of <50 in the mycelium and ≥500 in the fruiting body. According to the cutoff levels, 51 unigenes were selected as mature fruiting body-specific candidate genes; these unigenes were identified and shown to be involved in fruiting body maturation (Table 3). Overexpressed genes were selected as reference data based on read counts <50 in the fruiting body and ≥500 in the mycelium. At these cut-off levels, 172 unigenes were selected as mycelium-specific candidate genes that will be further examined to understand stage-specific expression of functional properties of L. edodes (Supplementary Table S9).
Among the 34 fruiting body maturation-specific DEGs, the LE_007730 unigene encodes an aspartic protease, which was expressed 609-fold (log 2 ratio = 9.23) higher in the fruiting body than in the mycelium (Table 3). Many proteases are used commercially, and fungal proteases have been reported in various mushrooms, including Agaricus bisporus, Armillariella mellaea, P. ostreatus, and Xylaria hypoxylon [49][50][51] . A novel aspartic protease from X. hypoxylon, a unique enzyme with relatively high thermostability, has been reported to inhibit HIV reverse transcriptase 51 . In addition, the gamma-glutamyl transpeptidase (LE_007464, log 2 ratio = 7.1), an important enzyme involved in flavor formation, is involved in the degradation of glutathione and glutathione-S-conjugates of xenobiotics and pharmaceuticals produced in detoxification processes 52 . Thus, gamma-glutamyl transpeptidase has been purified and characterized for further application in L. edodes 53 . Cyclohexanone monooxygenase (LE_007450, 8.34) is a valuable flavoenzyme that was initially isolated from Acinetobacter spp. and has been heterologously expressed for biocatalyst development 54 . Reverse transcriptase RNase H (LE_006898, 5.17) exhibits antifungal, antiviral, antitumor, and antiproliferative properties 8 . Ribonucleases, purified from the fruiting body of Russula delca, have been shown to exhibit antifungal and HIV-1 reverse transcriptase inhibitory activities. Aspartic protease (LE_007730), gamma-glutamyl transpeptidase (LE_007464), cyclohexanone monooxygenase (LE_007450), and reverse transcriptase RNase H (LE_006898) were expressed in the mature fruiting body stage only and not in the mycelium (Table 3). These results suggest that beneficial substances such as aspartic protease, gamma-glutamyl transpeptidase, cyclohexanone monooxygenase, and reverse transcriptase RNase H are produced optimally only during the mature fruiting body stage.
On the contrary, lectins are also abundantly expressed in the dikaryotic mycelium but not in the primordia of L. edodes 28 . Mushroom lectins have been reported to improve the antitumor and immunomodulatory activities of pharmaceutical substances 55 and to stimulate mushroom development 28,56 . Among the identified DEGs, the LE_008126 and LE_000173 unigenes, which encode lectin homologs, were upregulated 376-fold (log 2 ratio = 8.6) and 11.0-fold (3.5), respectively, in the mycelium compared with the fruiting body. The LE_008126 unigene in particular showed limited expression in the mature fruiting body (Supplementary Table S8). These results suggest that the expression of lectins is suppressed during fruiting body development as well as the primordia stage, and that lectins can be purified from mushroom mycelia to increase yield.
Among the identified 34 DEGs, the most significant expression difference was observed for the LE_007412 unigene, annotated as hydrophobin 1 (hyd1), which had 1,600-fold higher (log 2 ratio = 10.6) expression in the fruiting body than the mycelium ( Table 3). Expression of the hyd1 gene was limited in the dikaryotic mycelium, whereas the LE_000457 unigene, annotated as hydrophobin 2 (hyd2), was abundantly expressed in the dikaryotic mycelium in L. edodes, consistent with a previous report 28 . It has been reported that hydrophobins are essential for fruiting body formation in various mushrooms 28,57,58 , and hyd1 plays a more important role in fruiting body initiation (primordia) than maturation in L. edodes 57   only functions in fruiting body initiation but also in maturation in L. edodes. The results also indicate the validity of our transcriptome analysis. According to the stage-specific expression patterns of representative DEGs observed in this study, specific stages were verified as optimal for the production of beneficial nutritional and medicinal properties and new generation of improved L. edodes cultivars. Mycelium-specific DEGs were also listed in Supplementary Table S9. These useful substances are generally purified from mature L. edodes fruiting bodies and/or mycelia. No reports of comparative transcriptomic analyses of dikaryotic mycelia and mature fruiting bodies in L. edodes have been published. Thus, the large and accurate resource provided by our comparative transcriptomic analysis will contribute to increased production of L. edodes mushrooms with beneficial properties.
Finally, among 12 special DEGs for fruiting body maturation, DNA mismatch repair protein (LE_006999, 6.80), meiosis protein (LE_007462, 6.44), splicing factor (LE_011349, 6.59), ARM repeat-containing protein: RNA-binding protein (LE_001650, 3.92), heat shock protein 9 (LE_007440, 6.61), and cyclin: cell cycle control protein (LE_010885, 6.0) were involved in DNA replication, recombination, repair, and chromatin structure ( Supplementary Fig. S2 and Table 3). These proteins are involved in premeiotic replication, karyogamy, and meiosis, which occur during fruiting body maturation of L. edodes 18,47,48 . Interestingly, these results are consistent with the KOG and KEGG enrichment analyses of DEGs in this study. As shown in Supplementary Table S10, we explored 31 L. edodes genes related to functions including DNA replication, recombination, mismatch repair, spliceosomes, and chromatin structure and dynamics, which are important for premeiotic replication, karyogamy, and meiosis during fruiting body maturation. These results provide valuable genetic information for further studies of the molecular mechanism of fruiting body maturation in this species. The identification and analysis of DEGs of L. edodes provides valuable information regarding fruiting body maturation and the use of functional substances in future genetic studies.

Validation of transcriptome data by qRT-PCR.
Each DEG selected from L. edodes dikaryotic mycelia and mature fruiting bodies was analyzed by qRT-PCR. According to the KEGG annotation of genes differentially expressed between the dikaryotic mycelium and mature fruiting body, 36 DEGs from 18 KEGG subcategories in L. edodes were used to validate the results of the RNA-seq analysis (Fig. 8). One upregulated and one downregulated gene from each of the 18 KEGG subcategories were selected based on the significant fold changes in gene expression induced following formation of mycelia or mature fruiting body (Supplementary Table S11). The housekeeping gene 18S rRNA, used for qRT-PCR normalization 16 , was used as the reference gene, because the internal control gene was expressed at a constant level. These genes exhibited differential expression patterns similar to those in the RNA-seq experiments (Fig. 8). Therefore, DEGs identified by RNA-seq, which were validated by qRT-PCR, can be further investigated as candidate genes involved in the complex developmental stages of the mycelium and mature fruiting body.

Conclusions
This study is the first high-throughput transcriptome sequencing analysis of DEGs comparing the dikaryotic mycelium with the mature fruiting body of the commercially important Korean mushroom species L. edodes to investigate specific genes and specific developmental stages for optimal yield of beneficial properties and new generation. Furthermore, it is the first transcriptome analysis of functional enrichment using GO, KOG, and KEGG classification of DEGs to generate extensive and accurate genetic information on L. edodes for further study. Interestingly, functional annotations and enrichment analyses of the DEGs indicated that the genes overexpressed in the mature fruiting body stage are significantly involved in DNA replication, recombination, repair, chromatin structure, and the associated dynamics and significantly enriched in 'replication and repair' and 'transcription' pathways for premeiotic replication, karyogamy, or meiosis in L. edodes. Furthermore, many putative genes were verified by sequence comparisons with known the sequences of fungal gene families. Future analyses of the expression profiles of homologs related to fruiting body maturation and commercially useful functional substances would be useful to further explore the molecular mechanisms of L. edodes. DEGs showing significant changes in expression between different developmental stages can be further investigated for exploitation in commercial maturation and assessing the properties of the fruiting body using functional analysis and genetic engineering. Our transcriptome data can also be applied to transcriptome-based evaluation of the L. edodes stress response to mycoviral infection for research in an unexplored field.

Methods
L. edodes culture conditions. The commercial dikaryotic L. edodes strain, Sanjo701, was obtained from the Forest Mushroom Research Center, Korea. Mycelia were cultured from a source of L. edodes (mycelia discs, r = 5 mm) on potato dextrose agar at 25 °C in the dark. After 10 days, fresh mycelia covered the cellophane surface on top of the medium (diameter = 9 mm), and 10 g of mycelia were collected 59 . The mycelia were ground, rehydrated in culture broth, and inoculated inside the top of a sterile plastic sawdust bag with oak tree sawdust (850 g) and rice bran (150 g) with 55% moisture content. The sawdust bags with cotton aeration filters were incubated for 30 days at 25 °C in the dark to permit vegetative growth. To allow for mycelial ripening and browning, the incubated bags were transferred to a 20 °C greenhouse under natural light conditions for an additional 30 days. Next, the top piece of the bag, in which mycelia extended through more than half of the sawdust as a brown film, was cut and removed, and the bag was sprinkled with water for 12 h. The wet medium was maintained at 20 °C and >70% humidity to produce mature fruiting bodies on different sawdust bags. Fresh fruiting bodies were harvested immediately and stipes were removed to the extent possible 19 . Pilei with gills were frozen at −80 °C for use as representative samples of the mature fruiting bodies. L. edodes samples from the two developmental stages, the mycelium and fruiting body, were prepared for RNA extraction for further RNA-Seq and qPCR analysis.
SCIENtIFIC REPORTS | (2018) 8:8983 | DOI:10.1038/s41598-018-27318-z RNA isolation and cDNA library construction. Biological samples of total RNA were purified in triplicate from the mycelia and fruiting bodies of L. edodes as described previously 59 . RNA quality and quantity were evaluated using the Agilent 2100 BioAnalyzer with an RNA integrity value greater than seven (Agilent Technologies, CA, USA). We constructed a total of six libraries for RNA sequencing (RNA-Seq) using the Illumina TruSeq RNA sample preparation kit (Illumina, CA, USA), in accordance with the manufacturer's instructions. mRNA, which was enriched by poly A tail selection and chemically fragmented, was used for first-strand cDNA synthesis using a random hexamer primer, followed by second-strand cDNA synthesis. The cDNA fragments were blunted with adenine nucleotides and connected with sequencing adaptors. The libraries were then size-selected for the cDNA target fragments, and the selected cDNA was enhanced via polymerase chain reaction (PCR) using adaptor-specific primers. The cDNA library was quantified using the KAPA library quantification kit (Roche, Switzerland) according to the manufacturer's instructions.
Illumina sequencing and de novo assembly. High-throughput sequencing of each cDNA library from the mycelium and a mature fruiting body was performed in triplicate using the HiseqTM 2000 system (Illumina) to ensure the desired average sequencing depth. The 101-base pair (bp) raw paired-end reads were filtered according to the Phred quality score (Q ≥ 20) and read length (≥25 bp) using SolexaQA software (http://solexaqa.sourceforge.net/) to obtain high-quality clean reads by adaptor removal 30 . De novo assembly of the trimmed reads was performed using the Velvet (v1.2.10) 60 and Oases (v0.2.09) programs 61 , which implement the de Bruijn graph algorithm. We considered several hash lengths to determine optimal de novo assembly.
Raw sequencing data of three independent biological replicates for RNA-Seq were deposited in the National Center for Biotechnology Information (NCBI) sequencing read archive (https://submit.ncbi.nlm.nih.gov/subs/sra/) under accession number SRP137701. Figure 8. Validation of the RNA-seq results by quantitative real-time polymerase chain reaction. The log 2 fold change is denoted as the ratio of expression in the fruiting body to that in the mycelium. The y-axis represents the expression levels of the genes relative to the housekeeping gene (18S rRNA) for three biological and technical replicates. Significant differences of p < 0.01 were analyzed using Student's t-test.
SCIENtIFIC REPORTS | (2018) 8:8983 | DOI:10.1038/s41598-018-27318-z Short read mapping and expression profiles. Trimmed reads for each sequence tag were mapped to the assembled transcripts using Bowtie2 (v2.1.0) software 62 . The number of mapped reads for each unique transcript was calculated using an in-house script. Gene expression datasets were generated from each of two differential stages and raw counts were normalized and analyzed using the DESeq package in the R software 38 . The fold change (FC) and number of reads mapped to each unigene were used to identify DEGs between the two stages. A false discovery rate (FDR) was applied to calculate the p-value threshold for statistical significance in multiple-comparison tests. We used 'FDR ≤ 0.01 and | log 2 FC | ≥ 1' as the threshold to assess the significance of gene expression differences. Hierarchical clustering for all correlation analyses was performed using the Another Multidimensional Analysis Package library in R 39 . Upregulated and downregulated transcripts were subjected to Venn diagram analysis.
Functional annotation. All assembled transcripts from the total reads were validated by direct comparison using the fungi non-redundant (NR) (https://www.ncbi.nlm.nih.gov/refseq) and UniProtKB (http://www. ebi.ac.uk/uniprot) databases and the euKaryotic Orthologous Groups (KOG) tool (http://www.ncbi.nlm.nih. gov/KOG) using BLASTx (e-value ≤ 1e −10 ). Proteins with high sequence similarity were retrieved for analysis. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were carried out using the sequence similarities (e-value ≤ 1e −10 ) of proteins in the GO 33 and KEGG 63 databases. For each list of DEGs, down-and up-regulated genes were annotated using GO classification based on protein sequence similarity in the GO database. The number of genes assigned to GO terms was counted using in-house scripts from Seeders Co. (Korea).
Real-time quantitative (qRT) PCR analysis. cDNA was synthesized from approximately 1 μg of total RNA from the mycelium and mature fruiting body excluding the stipe using the GoScript Reverse Transcription system (Promega, WI, USA). SYBR green-based qRT-PCR was performed using the FastStart Essential DNA Green Master mix (Roche, Switzerland) using the LightCycler 96 system (Roche) to examine the expression levels of target and internal reference genes. Each qRT-PCR reaction was performed with 2 μl of first-strand cDNA and 100 nM each of forward and reverse primers in a total reaction volume of 20 μl. Detailed primer information is listed in Table S1. At least three independent biological and technical repeats were performed for each sample. The 18Sribosomal RNA (rRNA) gene was used as an internal control for normalization 16 . Relative gene expression levels were calculated by the comparative threshold cycle method using the LightCycler 96 qualitative detection module (Roche) 64 .