Abstract Lentinula edodes, one of the most highly regarded edible mushrooms in China, is susceptible to damage from high temperatures. However, a mutant strain derived from L. edodes, known as Le023M, has shown exceptional thermotolerance. Compared to the original strain Le023, Le023M exhibited accelerated mycelial recovery following heat stress. Through RNA-seq analysis, the majority of differentially expressed genes (DEGs) were found to be associated with functions such as “protein refolding”, “protein unfolding”, “protein folding”, and “response to heat”, all of which are closely linked to heat shock proteins. Furthermore, qRT-PCR results revealed significant accumulation of heat shock-related genes in Le023M under heat stress. GC-MS analysis indicated elevated levels of trehalose, aspartate, and glutamate in Le023M when subjected to heat stress. The highly expressed genes involved in these metabolic pathways were predominantly found in Le023M. Collectively, these findings highlight the following: (i) the crucial role of heat shock proteins (HSPs) in the thermo-resistant mechanisms of Le023M; (ii) the potential of trehalose accumulation in Le023M to enhance mycelium resistance to heat stress; and (iii) the induction of aspartate and glutamate accumulation in response to heat stress. These results shed light on the molecular mechanisms underlying the thermotolerance of Le023M, providing valuable insights for further understanding and improving heat stress response in L. edodes. The findings also highlight the potential applications of Le023M in the cultivation and production of L. edodes under high-temperature conditions. Keywords: Lentinula edodes, Heat stress, Protein folding, Metabolite, Trehalose 1. Introduction Lentinula edodes, which is known as the shiitake mushroom and is an edible mushroom in East Asia. It holds the distinction of being the major edible fungus, accounting for 22% of the world’s cultivated mushrooms [[37]1]. With its low fat content, L. edodes is rich in polyunsaturated fatty acids, constituting over 90% of its composition [[38]2]. Beyond its culinary appeal, L. edodes is highly regarded for its health benefits, thanks to active compounds like lentinan. Lentinan exhibits antimicrobial, anti-cancer, liver-protective, anti-diabetic, anti-flu, and anti-leukemia properties [[39][3], [40][4], [41][5], [42][6], [43][7]]. The agricultural cultivation and biotechnological applications of L. edodes have had a profound positive impact on people’s lives, owing to its valuable properties. In fact, China alone accounted for 90% of global production, yielding an estimated 11.64 million tons of L. edodes in 2020 ([44]https://www.chyxx.com/industry/202201/992805.html). The cultivation of L. edodes traces its roots back to China 800 years ago. Growth and development of edible fungi, including L. edodes, are significantly influenced by environmental factors, such as high temperatures. When L. edodes is cultivated at high temperatures, it experiences an increase in cell membrane fluidity, leading to structural changes in organelles and a decrease in enzyme activity [[45]8]. Mycelium development is negatively affected when temperatures exceed 32 °C, resulting in poor mycelia. Temperatures above 35 °C can even cause severe growth retardation or death. During the warmer months, the viability of spring-inoculated L. edodes mycelia deteriorates, rendering them vulnerable. Moreover, at high temperatures, the culture substrates or bags are more susceptible to pathogen infections, leading to reduced yields and compromised quality [[46]9,[47]10]. Based on transcriptomic, proteomic, and metabolomic analysis, the pattern of molecular expression of L. edodes in response to heat stress from mRNA, protein, and metabolites can be analyzed. This allows for a more comprehensive understanding of L. edodes molecular mechanisms in response to heat stress. It also provides insight into the development of strategies to improve its resistance to heat stress. The ion channel protein, G protein-coupled receptor, and histidine kinase are the first signal receptors to detect heat stress signals. Second messengers such as Ca^2+, ABA, and phosphoglycerol can then activate transcription factors. Then, these TFs activate the target genes involved in heat stress to produce a variety of protective proteins for cells, including heat stock proteins (HSPs) and late embryogenesis abundant (LEA) proteins [[48]11,[49]12]. HSPs serve as molecular chaperones and are synthesized in response to high temperatures, salinity, drought, starvation, and heavy metal ions [[50]13]. In abiotic stress, HSPs work as molecular chaperones, facilitating the refolding, stabilization, and assembly of other proteins without being part of a final structure [[51]14]. Trehalose is produced as a metabolite in response to heat stress. Heat stress induces endogenous trehalose accumulation in Pleurotus eryngii [[52]15]. Additionally, the supplementation of 100 or 200 g/L trehalose significantly inhibits the production of thiobarbituric acid-reactive substance (TBARS) in mycelial cells, providing cellular protection against superoxide toxicity during heat stress. One specific L. edodes strain, Le023M, was selected from the protoplasts of the original strain Le023, which had undergone UV treatment [[53]16]. Le023M was subjected to a heat shock at 37 °C for 24 h in potato dextrose (PDA) medium to evaluate its thermotolerance capacity and recovery rate. By conducting a comprehensive analysis of transcription and metabolism levels using integrated transcriptomic and metabolomic approaches, the mechanisms underlying heat stress were elucidated. This analysis provided valuable insights into the molecular expression patterns associated with heat stress responses. In response to heat stress, differentially expressed genes (DEGs) and differentially metabolites were identified and mapped into potential regulatory biological pathways. These pathways were crucial for the adaptation of the organism to its thermal environment. The findings from this study not only enhanced our understanding of the molecular foundations of heat stress but also paved the way for the discovery of previously unknown genes and pathways associated with heat response. This knowledge can be leveraged to develop strategies aimed at enhancing heat tolerance in crops and other organisms. 2. Materials and methods 2.1. Sample information Le023M, derived from the main cultivar of L. edodes in Hebei, Le023, is a mutant strain. Its authenticity was verified using an inter-simple sequence repeat (ISSR) marker, and it exhibited enhanced characteristics such as rapid growth and a pleasant taste [[54]17]. The culture technique findings demonstrated that Le023M exhibited a swift recovery from heat stress and promptly initiated the fruiting process. 2.2. Cultured and heat stress treatment Le023 and Le023M mycelia were taken out from 4 °C and activated three times on PDA medium. A 1 cm^2 section of mycelia-medium mixture was obtained by punching and transferred to a new PDA plate. The plate was then incubated for six days at 25 °C. Afterwards, the mycelia were cultured at 37 °C for 0 and 24 h ([55]Fig. 1A). In order to obtain optimal results, we conducted temperature screening and found that Le023 and Le023M exhibited the largest phenotypic difference at 37 °C for 24 h. Therefore, we selected 37 °C as the treatment temperature ([56]Fig. S1). Following that, all samples were transferred back to 25 °C and cultured for an additional 6 days. After the cultivation period, photos were taken of the samples. Fig. 1. [57]Fig. 1 [58]Open in a new tab Response to Heat Stress in Two Tested L. edodes Strains. (A) Culturing conditions for two strains. The mycelia were cultured on PDA medium at 25 °C for 6 days, followed by exposure to 37 °C for 0 or 24 h, and then cultured at 25 °C for an additional 6 days. (B) Phenotype of mycelia exposed to 37 °C for 0 and 24 h. (C) Observation of mycelia undergoing heat stress for 0 or 24 h. (D) Recovery of mycelial growth after heat stress in L. edodes. The values presented are the means ± standard deviation (SD) of three independent experiments. Statistical significance is denoted as **p < 0.01. Trehalose was added to the PDA medium until the appropriate final concentrations of 0, 5, 10, and 20 g/L were achieved. The medium was then sterilized at 121 °C for 20 min and cooled back down to 60 °C. Subsequently, the medium was separated into sterilized petri dishes. The heat treatment of mycelia was performed in the same manner. For microcopy analysis, a coverslip was inserted diagonally into a fresh PDA medium, positioned 3 cm away from the center of the Petri dish. Following incubation at 25 °C, allowing the mycelia to grow and spread onto the coverslip, the samples were subjected to a heat shock by transferring them to 37 °C for either 0 h (control) or 24 h. After the heat treatment, both Le023 and Le023M mycelia were incubated at 25 °C for additional six days, until the coverslip was fully colonized. Subsequently, a glass slide was placed on one side of the mycelia, which were then observed using a Leica STELLARIS 5 Confocal Microscopy (Leica Microsystems, Wetzlar, Germany). The resulting images were processed using Photoshop CS6 software. This approach allowed for detailed observation of the mycelial structures and collection of the required data. The activated mycelia were homogenized with medium using a homogenizer and then transferred into a triangular flask containing 100 mL of potato dextrose broth (PDB) medium. The cultivation process took place at 25 °C with agitation at a speed of 150 rpm for 7 days. Subsequently, 10 mL of the fungal endophyte culture supernatant was transferred to a new flask containing 100 mL of fresh PDB medium. This new mixture was further incubated under the same conditions, 25 °C and 150 rpm, for an additional 14 days. At the end of the 14-day period, the endophyte culture had reached a sufficient quantity for harvesting. The biomass was then separated from the broth and stored at −80 °C for further analysis. 2.3. Relative electric conductivity (REC) test Relative electric conductivity of the samples, abbreviated as REC, was measured as described in rice [[59]18]. Mycelia samples with different times of heat stress were placed into 50-mL tubes containing 30 mL of deionized water. The tubes were vibrated occasionally at 25 °C for 30 min, and then the electrical conductivity of the mixture was measured as E1. As a control, a blank consisting of only deionized water without sample (E0) was measured. After boiling the samples for 30 min, the electrical conductivity was measured as E2. The difference between E1 and E2 was used to determine the REC value. 2.4. Malondialdehyde (MDA) content The content of MDA was detected following the method of Heath and Packer [[60]19] with some modifications. Mycelia taken out from −80 °C were freeze-dried and ground into a fine powder. Triplicate 20 mg aliquots were extracted in 1 mL 10% of trichloroacetic acid (TCA) for 30 min at 4 °C. After centrifugation at 12,000g for 10 min, 0.2 mL of supernatant was transferred to a new 2 mL screwed centrifuge tube containing 0.8 mL of 0.67% of thiobarbituric acid (TBA). The mixture was boiled for 20 min. After centrifugation at 12,000g for 10 min, 300 μl aliquots in triplicates were transferred to 96-well flatted bottom plate. Absorbance was read at 532 nm and 600 nm, respectively. The content of MDA was calculated as follows: MDA content (nmol/g) = 53.763 × ΔA ÷ W (ΔA was difference of ΔA[532] and ΔA[600]. W was weight of mycelia). 2.5. H[2]O[2] content The hydrogen peroxide (H[2]O[2]) content was determined using the starch hydrogen peroxide assay kit (Megazyme) in a 96-well plate format. The assay was adapted from the manufacturer’s protocol, and appropriate dilutions were made. The results were expressed as the mean of three technical replicates, which were then displayed. 2.6. Measurement of ascorbate (AsA) and glutathione (GSH and GSSG) The AsA kit (Nanjing Jiancheng Bioengineering Institute, China) was used to determine the AsA content in accordance with the manufacturer’s instructions. The kit contains all the necessary reagents, including the AsA standard, for the AsA assay. GSH and GSSG were extracted from samples and the contents were measured using GSH kit and GSSG kit (Beyotime, China) according to the manufacturer’s instructions. The extraction process was carried out strictly following the guidelines provided by the manufacturer. 2.7. RNA extraction and library preparation Total RNA was extracted using the MagMAX™ mirVana™ isolation kit (Thermo Fisher), following the protocol of the manufacturer. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Samples with RNA Integrity Number (RIN) ≥ 7 were subjected to subsequent analysis. The libraries were constructed using TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. Then, these libraries were sequenced on the Illumina sequencing platform (HiSeqTM 2500), and 125 bp/150 bp paired-end reads were generated. 2.8. RNA-seq data analysis Raw data (raw reads) were processed using Trimmomatic [[61]20]. The reads containing ploy-N and low-quality reads were removed to obtain the clean reads. After removing the adaptor and low-quality sequences, the clean reads were assembled into expressed sequence tag clusters (contigs) and de novo assembled into transcript by using Trinity [[62]21] (version: trinityrnaseq_r20131110) in paired-end method. The longest transcript was chosen as an unigene based on the similarity and length of a sequence for subsequent analysis. The function of the unigenes was annotated by alignment with the NCBI non-redundant (NR), SwissProt, and clusters of orthologous groups for eukaryotic complete genomes (KOG) databases in Blastx [[63]22] with a threshold E-value of 10^−5. The proteins with highest hits to the unigenes were used to assign functional annotations. Based on the SwissProt annotation, GO classification was performed based on the mapping relationship between SwissProt and GO term. The unigenes were mapped to KEGG [[64]23] database to annotate their potential metabolic pathways. FPKM [[65]24] and read counts value of each unigene was calculated using bowtie2 [[66]25] and express [[67]26]. DEGs were identified using the DESeq package in Rstudio (version 4.2.1). Data were normalized using the function of estimate size Factors. P value and foldchange were determined using the function of nbinom Test. P value < 0.05 and foldChange > 2 or foldChange < 0.5 were set as the threshold for significant differential expression. Hierarchical cluster analysis of DEGs was performed to explore the transcript expression pattern. GO enrichment and KEGG pathway enrichment analysis of DEGs were performed using Rstudio (version 4.2.1) based on the hypergeometric distribution. Heatmap was generated using the pheatmap package in Rstudio (version 4.2.1). 2.9. Validation of RNA-seq data by qRT-PCR Quantitative reverse transcription-PCR (qRT-PCR) was performed to measure gene expression, following the procedure described by Mieog, Janeček [[68]27]. The TUB [[69]28] gene was selected as the reference gene for normalization. The target gene expression was determined by the E^−▵▵Ct method against the Ct value of TUB. The qRT-PCR primers used are listed in [70]Table S1. 2.10. Metabolome analysis Sample preparation, measurement, and analysis were performed as described by Zhao, Chen [[71]29]. The Agilent 7890B gas chromatography system, in conjunction with the Agilent 5977A MSD system (Agilent Technologies Inc., CA, USA), was used to analyze the derivative samples. To separate the derivatives, a DB-5MS fused-silica capillary column (30 m × 0.25 mm × 0.25 μm, Agilent J&W Scientific, Folsom, CA, USA) was utilized. A constant flow rate of 1 mL/min of helium (>99.999%) was employed as the carrier gas through the column. The injector temperature was maintained at 260 °C, and the injection volume was 1 μL in splitless mode. Initially, the oven temperature was set at 60 °C for 0.5 min. Subsequently, it was raised to 125 °C at a rate of 8 °C/min, then increased to 210 °C at a rate of 4 °C/min, followed by a further increase to 270 °C at a rate of 5 °C/min. Finally, the temperature was held at 305 °C for 3 min. The MS quadrupole and ion source (electron impact) were set to 150 °C and 230 °C, respectively. A collision energy of 70 eV was applied. Mass data was acquired in a full-scan mode within the m/z range of 50–500, with a solvent delay time of 5 min. To assess repeatability, quality control samples (QCs) were injected at regular intervals, specifically every 10 samples, throughout the analytical run. This provided a set of data for evaluating the consistency of the analysis. 2.11. Statistical analysis All experimental data were presented as mean ± standard deviation (SD) from a minimum of three replicates. Statistical significance between the datasets was evaluated using a t-test. GraphPad Prism 9.0 software was utilized for generating all figures. 3. Results 3.1. Heat stress inhibited growth of mycelia in L. edodes The effect of heat stress on L. edodes mycelia morphology was investigated by growing mycelia at 25 °C for 24 h. Heat stress substantially inhibited mycelia growth in both Le023 and Le023M ([72]Fig. 1B). Mycelia were dense and had more branches prior to heat stress ([73]Fig. 1C). Mycelia became sparse after heat stress and formed a visible heat stress circle. Microstructure study revealed that after exposure to heat stress, mycelia branches were diminished in both strains Le023 and Le023M ([74]Fig. 1C). Le023M mycelia grew faster than Le023Mycelia when cultivated at 25 °C, at 1.91 ± 0.28 mm/day for Le023 and 4.32 ± 0.19 mm/day for Le023M ([75]Fig. 1D). In two strains, heat stress inhibited mycelia growth. However, once heat stress was removed, Le023M mycelia resumed growth at a rate of 1.62 ± 0.22 mm/day for Le023M and 0.88 ± 0.11 mm/day for Le023. When mycelia were exposed to 39 °C for varying periods of time, Le023Mycelia died after 12 h, whereas Le023M did not cease developing until 24 h ([76]Fig. S2). This data suggested that both strains were impacted by heat stress, but mycelia of Le023M recovered quicker when heat stress was removed. 3.2. Effect on lipid peroxidation and non-enzymatic antioxidants The cell membrane may be damaged by high temperatures. Heat shock improved mycelia REC in both the Le023 and Le023M strains ([77]Fig. 2A). H[2]O[2] was accumulated in Le023M during heat stress, but there was a slight decrease in Le023 ([78]Fig. 2B). MDA content in Le023M was greater than in Le023 without heat stress. MDA level in Le023 increased in response to heat stress, but decreased in Le023M ([79]Fig. 2C). Heat stress caused a significant drop in AsA content in Le023. In contrast of Le023, Le023M showed a lower level of AsA with or without heat stress ([80]Fig. 2D). GSH level in Le023 increased swiftly in response to heat stress. Without heat stress, Le023M has a comparatively high amount of GSH compared to Le023. GSH content decreased dramatically after heat stress ([81]Fig. 2E). Le023 showed an extremely high GSSH concentration prior to heat stress. Le023 and Le023M demonstrated a drop in GSSG content after heat stress and maintained a very low level ([82]Fig. 2F). Fig. 2. [83]Fig. 2 [84]Open in a new tab Physiological parameters of mycelia during heat stress. (A) Relative electrolytic conductivity (REC) of mycelia under heat stress. The content of superoxide (B), MDA (C), AsA (D), GHS (E), and GSSH (F) in mycelia under heat stress. Three biological replicates were conducted for each measurement. Asterisks indicate significant differences between the two strains (**p < 0.01). 3.3. Differentially expressed genes (DEGs) under heat stress The differentially expressed gene analysis was conducted to investigate the response to heat stress in the mycelia of Le023 and Le023M strains using transcriptome sequencing on the Illumina Hiseq-2000 system. The transcriptome profiles of Le023 and Le023M strains' mycelia were analyzed after subjecting them to 37 °C for 0 and 24 h. Following the removal of adapter sequences and low-quality reads, an average of 48,070,608 clean reads and 7,031,504,121 clean bases were obtained ([85]Table 1). The average percentage of valid bases and Q30 scores (Phred score ≥ 30) were 95.96% and 94.46%, respectively. By aligning the clean data without a reference genome, a total of 16,139 unigenes with a combined length of 28,076,978 bp were identified. Table 1. RNA-seq data summary and de novo transcriptome assembly statistics for Le023 and Le023M. Sample clean_reads clean_bases valid_bases Q30 GC Le023 0 h 48057706 7034013881 94.93% 94.36% 48.95% Le023 24 h 47988430 7023711031 94.94% 94.37% 48.84% Le023M 0 h 48320358 7074662299 95.10% 94.56% 49.06% Le023M 24 h 47915940 6993629271 94.88% 94.52% 48.81% [86]Open in a new tab Comparing the gene expression levels to the control (0 h of heat stress), a total of 1421 unigenes in the Le023 strain showed significant differential expression (more than two-fold) under heat stress ([87]Fig. 3A and E). Among these DEGs, 833 were upregulated, while 588 were downregulated in response to heat stress ([88]Fig. 3E). Similarly, in the Le023M strain, a total of 1483 DEGs were identified under heat stress ([89]Fig. 3B and E), with 517 unigenes upregulated and 966 unigenes downregulated ([90]Fig. 3E). Importantly, significant differential expression was observed between the Le023 and Le023M strains, with 1495 and 1643 unigenes showing differential expression in the absence or presence of heat stress, respectively ([91]Fig. 3C–E). Following the heat stress, a total of 231 heat-induced DEGs and 277 heat-repressed DEGs were observed in both Le023 and Le023M strains ([92]Fig. 3F and G). These results suggest that the Le023 and Le023M strains display distinct responses to heat shock. Fig. 3. [93]Fig. 3 [94]Open in a new tab Analysis of Differentially expressed genes (DEGs) between two strains in response to heat stress. Volcano plots of DEGs in different comparisons: (A) Le023 24 h versus Le023 0 h; (B) Le023M 24 h versus Le023M 0 h; (C) Le023M 0 h versus Le023 0 h; (D) Le023M 24 h versus Le023M 24 h. In the volcano plots, red, green, and gray dots represent upregulated genes, downregulated genes, and genes without significant difference, respectively. (E) Column diagram of DEGs in different comparisons. Pink and azure column represent upregulated and downregulated genes, respectively. Y-axis indicates the number of genes. Venn graph for Le023 and Le023M based on heat-induced genes (F) and heat-repressed genes (G). (For interpretation of the references to