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).