Comparative transcriptome analysis on candidate genes associated with fruiting body growth and development in Lyophyllum decastes

Sample preparation

The L. decastes strain was preserved at the China General Microbiological Culture Collection Center. The mycelia were cultured on potato dextrose agar with a pH of 4–5 at 25 °C in the dark. The cultivation substrates, contained in plastic bags, were comprised of 750 g medium consisting of 75% cottonseed hull, 12% sorghum seed, 12% wheat bran, 0.5% sugar, and 0.5% superphosphate. The bags were incubated in a 19 °C greenhouse under artificial light conditions once fully colonized. L. decastes samples were collected from randomly selected bags at five different developmental stages: the vegetative mycelium stage (VM), primordium initiation stage (PI), young fruiting body stage (YF), medium-size fruiting body stage (MSF), and mature fruiting body stage (MF) (Fig. 1). The morphological changes were observed, which the VM turned into PI in 7–10 days, and following young fruiting body appeared (YF) in 3–5 days, and the medium-size fruiting body appeared (MSF) in 5–7 days, and finally fruiting body matured (MF) in 5–7 days. All samples were promptly cryopreserved in liquid nitrogen and stored at −80 °C for RNA extraction.

RNA extraction and sequencing

The acquisition of total RNA from each sample was executed through the utilization of diverse kits, including the TRIzol Reagent (Invitrogen, Waltham, MA, USA). After extracting the RNA, its quality and quantity were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), NanoDrop (Thermo Fisher Scientific Inc. Waltham, Ma, USA), and 1% agarose gel. For each sample, 1 μg of total RNA with a RNA integrity (RIN) value greater than 7 was utilized for subsequent library preparation. The NEBNext First Strand Synthesis Reaction Buffer and NEBNext Random Primers were used for mRNA fragmentation and priming. ProtoScript II Reverse Transcriptase was utilized to produce the initial strand cDNA, followed by the Second Strand Synthesis Enzyme Mix to generate the complementary second strand cDNA. The size selection of Adaptor-ligated DNA was carried out using AxyPrep Mag PCR Clean-up (Axygen, Union City, CA, USA) and fragments measuring approximately 360 bp (with an approximate insert size of 300 bp) were recovered. The manufacturer’s protocol (NEBNext^® Ultra^™ RNA Library Prep Kit for Illumina^®) was followed in order to construct next-generation sequencing library preparations. The Illumina HiSeq machine was used according to the manufacturer’s instructions (Illumina, San Diego, CA, USA) to load and multiplex the libraries with distinct indices. The statistical power of this experimental design calculated in RNASeqPower was 0.87. Three biological replicates were taken for each stage.

Assembly and functional annotation

The initial paired-end reads underwent a series of quality control measures, which involved the elimination of adapter, low-quality, and duplicate sequences. The quality of the reads was evaluated by FastQC software (Version 0.10.1). In order to remove technical sequences, including adapters, polymerase chain reaction (PCR) primers, or fragments thereof, and quality of bases lower than 20, pass filter data of fastq format were processed by Cutadapt (Version 1.9.1, phred cutoff: 20, error rate: 0.1, adapter overlap: 1 bp, min. length: 75, proportion of N: 0.1) to be high quality clean data. Subsequently, a de novo assembly was conducted using Trinity (version 2.4.0) (Grabherr et al., 2011), a novel approach that facilitated the efficient and robust reconstruction of transcriptomes from RNA-seq data. This approach employed three independent software modules, namely Inchworm, Chrysalis, and Butterfly, in a sequential manner to process large quantities of RNA-seq reads. The contigs obtained were further refined using Cd-hit (version 4.6.2) to eliminate duplicates, resulting in the production of the unigene sequence file. To annotate the assembled unigene sequences, they were searched against several databases, including non-redundant protein of NCBI (Nr, https://www.ncbi.nlm.nih.gov/protein), Gene Ontology (GO, http://www.geneontology.org), Cluster of Orthologous Groups (COG, http://www.ncbi.nlm.nih.gov/COG/), Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp/), and Swiss-Prot (http://www.expasy.org), using blast software (version 2.2.28+) with an E-value of less than <10⁻⁵.

Expression and enrichment analysis

The reference gene file utilized in this study was the unigene sequence file. Using RSEM (v1.2.6), gene and isoform expression levels were estimated from the pair-end clean data, according to the methodology proposed by Fu et al. (2012). In order to identify the genes that were differentially expressed, the DESeq2 (version 1.26.0) Bioconductor package was utilized (Anders, McCarthy & Chen, 2013; Love, Huber & Anders, 2014), which employed a model based on the negative binomial distribution. The adjustment for controlling the false discovery rate (FDR) was carried out using the Benjamini and Hochberg approach, the detection of DEGs was set at a significance level of FDR < 0.05 and |log₂ fold change| > 1. The genes that were enriched were determined using GO-TermFinder (version 0.86), a tool that assigned annotations to a set of genes based on a statistically significant p-value less or equal than 0.05. Additionally, we developed in-house scripts to identify significant DEGs in KEGG pathways. Pathway significance enrichment analysis using KEGG pathway as unit, hypergeometric test was applied to find out the pathway of significant enrichment in DEGs compared with the reference background.

qPCR

To ensure the accuracy of the RNA-Seq data, a representative selection of DEGs was chosen. These genes were involved in the development of L. decastes fruiting bodies and encoded cellulose, galactosidase, and laccase. The fragments per kilobase of exon model mapped fragments (FPKM) values for these genes are available (Table S8), while the primer sequences are listed in Table S9. RNA extraction was carried out using RNAios Plus manufactured by Takara, and cDNA was synthesized using the Revert Aid First Strand cDNA Synthesis Kit provided by Thermo Scientific. The transcriptome data was utilized to obtain specific primer pairs, which were subsequently employed for qRT-PCR. The CFX96 TouchTM Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) was used to conduct the qRT-PCR in triplicate. The reaction mixture consisted of 20 μL, which included 10 μL of 2×SYBR Green mixture (TaKaRa, Dalian, China), 1 μL of diluted cDNA, and 1 μL each of specific forward and reverse primer. The 2^−ΔΔCT method was used to analyze the relative gene expression, with β-Actin serving as an internal control (Deepak, Ishii & Park, 2008).

Determination of extracellular enzyme activity

The culture medium was sampled at the VM, PI, YF, MSF and MF stages of L. decastes in order to determine extracellular enzyme activity. The activities of cellulase, galactosidase, laccase and peroxidase were determined by Lanzhou Boruke Biotechnology Co., LTD., China.

Illumina sequencing and assembly

The objective of this study was to identify changes in the transcriptome of L. decastes at five distinct developmental stages. To achieve this, RNA samples were obtained from the VM, PI, YF, MSF, and MF stages, which were then subjected to high-throughput sequencing. Overall, approximately 63 Gbp of raw data were produced from 15 samples across the VM, PI, YF, MSF and MF stages (Table S1). After data filtering and trimming, a total of 743,342,738 high-quality clean reads were generated with a Q20 percentage >97.3% (Table S2). Utilizing Trinity software, the reads were assembled into 54,338,739 contigs with a mean length of 86.2 bp and an N50 length of 103 bp (Table S3). The contigs were initially assembled, and a further assembly process was subsequently conducted, resulting in the formation of 15,451 unigenes. According to a recent research, 14,499 gene models have been predicted from the genome of L. decastes (Xu et al., 2023). Compared with these previously predicted genes, the average GC content of the unigenes in this research (51.78%) is simililar to that of L. decastes LGR-d1-1 (50.25%) and LRG-d1-5 (50.38%). Moreover, we de novo assembled some putative novel unigenes, and will provide important information for further annotation and characterization of L. decastes. The length of these unigenes ranged from 201 bp to 24,176 bp, with an average length of 1,462.68 bp and N50 length of 2,921 bp (Table S3). Of all the unigenes, 6,679 (43.23%) were between 200 to 1,000 bp long, 6,280 (40.64%) were within the range of 1,000~3,000 bp in length, while 16.13% were longer than 3,000 bp. Notably, none of the unigenes were shorter than 200 bp (Fig. 2). In order to enhance the comprehensiveness of the L. decastes transcriptome sequencing, the assembled unigenes were utilized to map the clean reads. The mapping outcomes are presented in Table S2.

Sequence annotation

The compiled unigenes were subjected to annotation through the utilization of the BLASTx algorithm against sequences sourced from four publicly available gene databases. A total of 11,783 (76.84%) unigenes were annotated in at least one database, 10,907 (70.59%) unigenes were observed to matched known proteins in the Nr database (E-value < 1.0E⁻⁴⁵), followed by the KEGG (6,650, 43.04%), SwissProt (6,454, 41.77%) and COG (5,324, 34.46%) database (Table S4 and Fig. S1). However, only 912 (5.90%) unigenes matched in all database. After functional annotation, the species distribution analyses of the annotated unigenes were calculated from the annotation results. We found that the top ten representative annotated species in Nr database were L. decastes LRG-d1-1 (66.13%), H. marmoreus (25.45%), Termitomyces sp. J132 (7.98%), Galerina marginata CBS 339.88 (1.93%), Laccaria amethystine LaAM-08-1 (0.95%), Pelomonas sp. Root 1217 (0.84%), Hebeloma cylindrosporum h7 (0.78%), Laccaria bicolor S238N-H82 (0.69%), Piloderma croceum F1598 (0.65%) and Pelomonas sp. Root1444 (0.63%), as shown in Fig. 3. Significantly, BLAST comparison was performed between the 15,451 unigenes in this study and 14,499 genes in the study by Xu et al. (2023), which revealed that 74.25% unigenes can be found in the 14,499 genes with an identity of 85%.

In order to evaluate the comprehensiveness of the transcriptome and the effectiveness of the annotation process, an investigation was carried out on the annotated unigenes involved in GO classifications. The three primary GO categories, namely biological process, cellular component, and molecular function (Fig. S2), made up a total of 4,082 unigenes that were classified into 59 functional groups. The matched unigenes were allocated to 25 classifications by the COG database, with the most extensive category being exclusively related to general function prediction and encompassing 5,324 unigenes (Fig. S3). Moreover, the unigenes that underwent annotation in the KEGG database were categorized into six distinct biological pathways, which include cellular process, environmental information processing, genetic information processing, human diseases, metabolism, and organismal system (Fig. S4).

Analysis of differentially expressed genes (DEGs)

To identify and evaluate the transcriptional profiles during L. decastes fruiting body formation and development, the expression levels of unigenes from five developmental stages were calculated using the Bioconductor software package and DESeq2 (Anders, McCarthy & Chen, 2013; Love, Huber & Anders, 2014), setting a FDR value ≤ 0.05 and log₂ ratio ≥ 1. Pairwise comparisons showed that 2,502, 713, 208 and 1,449 unigenes were significantly differentially expressed between PI and VM, YF and PI, MSF and YF, and MF and MSF, respectively, and 1,568, 1,673 and 2,704 unigenes were significantly differentially expressed between MF and YF, MF and PI, and MF and VM, respectively (Figs. 4, S5, and Table S5). As shown in Fig. 4, most DEGs occurred between the MF and VM (1,278 up-regulated, 1,426 down-regulated) and PI and VM (1,185 up-regulated, 1,317 down-regulated), and only a minimal number of unigenes were identified as DEGs between MSF and YF (70 up-regulated, 138 down-regulated), illustrating that a significant number of genes were responsible for the conversion of mycelium into fruiting body. When compared to the initial phase of fruiting body development, a higher quantity of genes displaying significant alterations was detected during the subsequent stage. This suggested that the maturation period of the fruiting body is a complex and multifarious process. As predicted, there was a gradual increase in the number of differentially expressed genes (DEGs) during the transition from mycelium to a fully matured fruiting body. Notably, in all compared groups, the quantity of unigenes that were down-regulated surpassed the quantity of those that were up-regulated. To investigate gene expression relationships, Venn diagrams of different pairwise comparisons were introduced to show the statistics of the DEGs in depth (Fig. 5). In total, only 19 unigenes were differentially expressed both in PI vs VM, YF vs PI, MSF vs YF and MF vs MSF. However, 323 DEGs were detected in MF vs MSF, MF vs YF, MF vs PI, and MF vs VM, simultaneously. These genes, which were significantly differentially expressed at each stage, may play a crucial role in the development of the L. decastes fruiting body.

Some main up-expressed and down-expressed DEGs were detected between different development stages of L. decastes (Table 1). We found that the expression difference of these genes in PI vs VM were significantly higher than those in other comparison groups. Of these unigenes, with log₂ fold change values ranging from −9.42 to 10.06 some were regarded as new genes that encoded unknown proteins or hypothetical proteins that could not be matched in any database. In different comparisons, some unigenes showed different up-regulation or down-regulation and need to be further studied.

Gene ontology (GO) enrichment analysis of DEGs

With the objective of identifying the potential contributions of genes that were differentially expressed, we carried out a GO functional enrichment analysis to assign their functions to the process of fruiting body development. Subsequently, the DEGs were subjected to GO term enrichment analysis, which enabled the characterization of the relevant GO terms. The results showed that the DEGs exhibited connections with a multitude of biologically significant activities. The DEGs were sorted into three main GO categories, with each consisting of 30 distinct functional subcategories. The number of DEGs in each subcategory was determined using the GO database, and the count ranged from 1 to 329 with Log₁₀ p-vlaue > 1. The significance threshold for the analysis was set at a corrected p-value of ≤0.05. For different comparison groups, the DEGs were significantly enriched in 26, 23, 20, 23, 25, 26 and 27 GO categories in PI vs VM, YF vs PI, MSF vs YF, MF vs MSF, MF vs YF, MF vs PI, and MF vs VM, respectively (Figs. 6 and S6). After analyzing the seven comparison groups, a significant enrichment of the categories “metabolic process” (GO:0008152), “cell part” (GO:0044464), and “catalytic activity” (GO:0003824) within the domains of biological process, cellular component, and molecular function, respectively, was observed. Furthermore, a substantial proportion of DEGs were linked to functional categories such as “cellular process” (GO:0009987), “single-organism process” (GO:0044699), “membrane part” (GO:0044425), and “binding” (GO:0000988). A certain number of DEGs belonged to the biological process category “response to stimulus” (GO:0050896) and molecular function category “antioxidant activity” (GO:0016209) in different comparison groups except for comparison MSF vs YF comparison. Analogously, no DEGs in MSF vs YF and YF vs PI were annotated to molecular function categories “developmental process” (GO:0003006), multi-organism process (GO:0051704) and “immune system process” (GO:0002376).

Our research focused on investigating the molecular mechanisms involved in the formation and development of fruiting bodies. To accomplish this, we conducted a GO functional classification analysis of DEGs at various stages of development. Our analysis revealed that the up-regulated DEGs in the PI vs VM and MF vs VM comparisons were primarily associated with the GO categories of “cellular process” (GO:0009987), “single-organism process” (GO:0044699), “metabolic process” (GO:0008152), and “binding” (GO:0000988). Conversely, there were relatively few down-regulated DEGs in these comparisons. However, in other comparisons, the opposite trend was observed (Figs. 6 and S6). These findings suggest that the expression profiles of genes associated with these GO categories may differ between the mycelium and fruiting body stages of L. decastes development.

Pathway enrichment analysis of DEGs

To further scrutinize the biological functions of the DEGs and synthesize the biological pathways that pertain to the fruiting process, we mapped the DEGs from each comparison to the KEGG database (Figs. 7, S7 and Table S7). Interestingly, in the five developmental stages of L. decastes, the DEGs were only enriched in five pathways in YF vs PI, and there were 13, 17, and 16 pathways in PI vs VM, MSF vs YF and MF vs MSF, respectively. The KEGG sub-classification analysis revealed that a significant number of DEGs were implicated in metabolic pathways (ko01100), ribosome (ko03010), spliceosome (ko03040), biosynthesis of second metabolites (ko01110) and RNA transport (ko03013). These pathways are associated with genetic information processing and metabolism, and their prominence in the comparison between PI and VM indicates that they may play a pivotal role in the conversion of mycelium to primordium. Nevertheless, we found no DEGs enriched in the pathway of glycerolipid metabolism (ko00561) in the comparison of PI vs VM, but found them in the rest of the comparisons. In addition, there were 16, 15 and 13 DEGs assigned to the pathway glycosphingolipid biosynthesis-globo (ko00603) in PI vs VM, YF vs PI and MSF vs YF, respectively, but not in MF vs MSF. The comparative analysis of MSF vs YF, MF vs MSF, MF vs YF and MF vs, VM revealed the presence of DEGs that were enriched in several pathways related to environmental information processing. These pathways included the MAPK signaling pathway (ko04013), FoxO signaling pathway (ko04068), HIF-1 signaling pathway (ko04066), TNF signaling pathway (ko04668), and VEGF signaling pathway (ko04370). These results suggest that various signaling pathways may play crucial roles in the transition from a young fruiting body to mature fruiting body.

Meanwhile, we focused on the specific DEGs that were enriched in different pathways across multiple comparisons. We found that TRINITY_DN24453_c5_g1_i2, a cell cycle checkpoint control protein galactosidase, was simultaneously enriched to 14 metabolic pathways of PI vs VM, YF vs PI, MSF vs YF and MF vs MSF. Similarly, three other DEGs, TRINITY_DN21188_c3_g1 (hypothetical protein), TRINITY_DN23599_c5_g3 (polyketide synthase) and TRINITY_DN24271_c1_g2 (hypothetical protein) were simultaneously identified in the various pathways of MF vs MSF, MF vs YF, MF vs PI and MF vs VM.

Special genes involved in fruiting body development

Previous investigations found that the analysis of DEGs offered valuable information regarding the development of the fruiting body. This process is initiated by a stimulus that activates particular genes responsible for the transformation of mycelium into a fruiting body (Yu et al., 2012; Muraguchi et al., 2015; Song, Kim & Kim, 2018; Yan et al., 2023). To identify the specific unigenes that contribute to the development of fruiting bodies in L. decastes, we conducted an analysis of the gene expression profiles across different comparisons. Our results indicated that 41 extracellular enzyme genes exhibited significant differences in expression levels at various stages, including eight cellulase genes, 11 galactosidase genes, 14 laccase genes, and eight peroxidase genes (as outlined in Table 2). Among these DEGs, there were 29, 14, three and 24 in PI vs VM, YF vs PI, MSF vs YF and MF vs MSF, respectively. The statistics also showed that a total of five, six and six extracellular enzyme genes were exclusively differentially expressed in PI vs VM, YF vs PI and MF vs MSF, respectively. In particular, the unigene TRINITY_DN24453_c5_g1_i2 that encoded a galactosidase was found to be significantly down-regulated in PI vs VM and MSF vs YF and up-regulated in YF vs PI and MF vs MSF. However, none of the laccase and peroxidase genes were differentially expressed in MSF vs YF or MF vs MSF.

Additionally, our primary emphasis was directed toward DEGs that fell under the category of transcription factor families. As a result, we systematically compiled a list of 19 potential transcription factors that were categorized into 14 distinct types. All of these types exhibit noteworthy variations in their expression levels throughout diverse developmental stages of L. decastes, as demonstrated in Table 3. Among of them, different amounts of several types transcription factors were characterized, including bHLH, C2H2-type, MADS box, MYB, and TFIIA, etc. While significant changes in expression were observed in these transcription factor genes, their patterns of up-regulation or down-regulation varied across different comparisons. For instance, the MYB gene TRINITY_DN24473_c6_g6_i2, was down-regulated in PI vs VM and MSF vs YF, up-regulated in YF vs PI and not detected in MF vs MSF. In summary, the transcription factor genes manifested dynamic expression profiles while the fruiting body of L. decastes was forming.

Our analysis led us to conclude that the extracellular enzyme genes and transcription factor genes found in L. decastes played a role in the initiation and maturation of fruiting bodies. Moreover, the discovery and examination of the distinctive genes of L. decastes provide significant genetic knowledge for future investigations into the molecular processes involved in the growth of fruiting bodies.

Validation of RNA-Seq data by qRT-PCR

In order to validate the transcriptome data, we analyzed specific genes encoding extracellular enzymes that were associated with the development of fruiting bodies using qRT-PCR to determine their expression patterns. The results were in agreement with those obtained by RNA-Seq. The conformity between the two approaches affirmed the trustworthiness of our RNA-Seq data (Fig. 8).

Analysis of extracellular enzyme activity

By releasing extracellular enzymes, fungi can effectively transform complex macromolecules into smaller units, thereby providing the necessary nutrients for their growth and development. Through the above transcriptome data analysis, we found that many extracellular enzyme-related genes showed significant differential expression during the fruiting body development of L. decastes (Table 2 and Fig. 8). To investigate the alterations in extracellular enzyme activity of L. decastes during various stages of growth and development, we established a theoretical framework to comprehend the utilization mechanism of L. decastes in a medium. We analyzed the activities of cellulase, galactosidase, laccase and peroxidase across five developmental stages (VM, PI, YF, MSF, and MF) of L. decastes. The experimental results showed that all four extracellular enzymes showed high activity in primordial initiation stage, and then decreased gradually (Fig. 9). The activities of cellulase, galactosidase and laccase all reached the maximum values at the primordium initiation stage, while the activities of peroxidase reached the maximum value at the vegetative mycelium stage (Fig. 9). The enzyme activity experiments conducted on L. decastes demonstrated that the extracellular enzymes exhibited discernible differences across various stages of growth and development. In order to elucidate and confirm the differences of these extracellular enzymes in different stages, expression levels of three cellulase regulators were selected and analyzed by qRT-PCR (Fig. S8). We found that the expression level of transcript inhibitors CRE1 and ACEI was significantly down-regulated at the primordium initiation stage, and then obviously up-regulated at the young fruiting body stage. The expression difference in the patterns of transcription activator ACEII is similar to that of cellulase activity between different developmental stages. The subsequent examination of enzyme activity provided essential insights into the pivotal role that extracellular enzymes play in the growth and development of the fruiting body in L. decastes.

Extracellular enzyme-related genes are important for the formation of L. decastes fruiting body

In edible fungi, the substrates used for cultivation are usually composed of wheat bran, cottonseed husk, sawdust, and other substances that possess significant nutrient molecules, including lignin, cellulose, hemicelluloses, starch, and protein (Huang et al., 2019). To meet the demands of growth and metabolism, edible fungi are stimulated to produce hemicellulose, cellulose, and various extracellular hydrolases during their development. These enzymes are essential for the degradation, assimilation, conversion, and utilization of large nutrient molecules (Lechner & Papinutti, 2006; Rani, Kalyani & Prathiba, 2008; Barh et al., 2022; Guigón-López et al., 2014). For example, Volvariella volvacea produces multiple forms of extracellular laccase when grown in specific culture conditions, and some laccase genes may play an important role in the development of fruiting bodies (Chen et al., 2003; Lu et al., 2015). Studies have found that amplifying the operation of laccase stimulates the initiation of H. marmoreus primordia and elevates the production of fruiting bodies (Zhang et al., 2015b; Zhang et al., 2016). Organisms’ carbohydrate metabolismrely heavily on lignocellulolytic enzymes, which are crucial for the transformation of the lignocellulolytic matrix during fungal growth and serve as a necessary physiological function (Elisashvili et al., 2002; Bánfi et al., 2015). The peak of extracellular cellulase activity was positively correlated with the biomass of fruiting bodies in Agaricus bisporus (Claydon, Allan & Wood, 1988). Moreover, the importance of peroxidase in the growth of mycelium and the development of fruiting bodies in H. marmoreus, Botrytis cinerea, and Ganoderma lucidum has been demonstrated (Kim et al., 2011; Ren et al., 2017; Chen et al., 2018).

In the current study, a group of genes encoding extracellular enzymes such as cellulose, galactosidase, laccase and peroxidase were up-regulated or down-regulated from mycelium to mature fruiting body stages (Table 2). The analysis of qPCR and enzymeactivityalso indicated that the gene expression and enzyme activity of extracellular enzymes showed significant differences at different developmental stages (Figs. 8 and 9). Surprisingly, the differences in enzymatic activity and the expression levels showed different patterns of change, such as galactosidase activity and galactosidase genes (Figs. 8 and 9). We hypothesize that the process of extracellular enzymes affecting fruiting body development is controlled by a complex regulatory network involving multiple genes, and futher work is need to characterize the mechanism of these extracellular enzymes and related genes in fruiting body development. Specifically, the gene encoding galactosidase (TRINITY_DN24453_c5_g1_i2) was not only significantly differentially expressed in each comparison but was categorized into 14 KEGG pathways. The outcomes of this investigation are consistent with earlier studies, indicating that the aforementioned genes are likely linked to the regulation of fruiting body formation in L. decastes through direct or indirect involvement, but further research is needed. Therefore, additional investigation is essential in order to recognize and evaluate the existence and roles of these extracellular enzyme genes in the growth of L. decastes mushrooms.

Transcription factors are essential for the formation of L. decastes fruiting body

Transcription factors act as molecular switches, binding to specific DNA sequences to regulate gene expression, and participating in many important cellular processes. In fungi, transcription factors that regulate fruiting body formation have been identified (Ohm et al., 2011; Park, 2013; Xiang et al., 2014; Hamann, Osiewacz & Teichert, 2022). In this research, we identified some different types of transcription factors that might play important roles in fruiting body formation, such as C2H2 type zinc finger, GATA zinc finger, HMG-box, MADS-box and TFIIA (Table 3). In filamentous fungi, C2H2 zinc finger protein is critical for the regulation of carbon metabolism and the expression of cellulase and hemicelluloses which are related to cellulose degradation, such as Trichoderma reesei cellulase regulators CRE1, ACE, and Neurospora crassa 11 hemicellulose regulators XLR-1(Martinez et al., 2008; Sun et al., 2012; Campos Antonieto et al., 2014). It has been found that exp1, a HMG-box transcription factor is a crucial factor in fruiting body morphogenesis and takes part in the regulation of fruiting body opening and autolysis (Muraguchi et al., 2008). Previous research has determined that there are two L. edodes transcription factors, PriA and PriB, that are involved in the modulation of gene expression during the process of primordium formation (Endo et al., 1994; Miyazaki et al., 2004). In addition, during the development of fruiting bodies in S. commune, there was an up-regulation of fst4, a transcription factor that is specific to fungi (Ohm et al., 2011). Recently some transcription factors including MADS, C2H2, and FST4 were identified in the mushroom formation of P. tuoliensis (Fu et al., 2017). Therefore, it is postulated that the transcriptional regulators coded by these unigenes are noteworthy genetic resources for forthcoming studies on the inception and maturation of primordia, along with the formation of fruiting bodies in L. decastes. However, it needs more data and experiment areas to prove it.

Signaling pathways in L. decastes fruiting body development

KEGG pathway analysis is a useful way to explore the metabolic and regulatory networks in organisms and is often used to search for pathways with significant enrichment in DEGs. According to earlier investigations, a significant number of KEGG pathways were pertinent to the formation of fruiting bodies in fungi, which included primary carbohydrate metabolism and the MAPK signaling pathway (Ceccaroli et al., 2011; Fu et al., 2017). The KEGG pathway analysis of different developmental stages of L. decastes revealed that a significant number of DEGs were enriched in various signaling pathways, as illustrated in Figs. 7 and S7. This observation is consistent with previous studies on other species, highlighting the correlation between signaling pathways and fungal fruiting body development. For instance, IDC2 and IDC3 genes were involved in signaling pathways and linked to Podospora anserine fruiting body development (Lalucque et al., 2017). Similarly, the MAPK signaling pathway was found to be up-regulated during fruiting body development in Coprinopsis cinerea (Cheng et al., 2013). However, studies about MAPK signaling L. decastes in were lacking. It should be emphasized that the DEGs involved in signaling pathways were exclusively associated with processing environmental information. This implies that environmental cues may significantly impact the development of fruiting body morphology. These findings collectively indicate that signaling pathways are essential for L. decastes fruiting body formation. The elucidation of the molecular mechanisms involved in L. decastes fruiting body development necessitates further investigation into the DEGs mentioned above. Consider that the current work involved only the transcriptome data, it is only possible to predict some possible functions of identified DEGs, and their specific regulatory mechanisms and networks need to be further studied and clarified.