The two sexually compatible monokaryotic strains SP3 and SP30 were obtained from the dikaryotic L. edodes strain JZB2102217, using protoplasting method as described by Gao et al. (2022). The strain JZB2102217 is a widely used high-temperature type cultivar in northern China. Mycelia from original cultures were punched out using a cork borer (1 cm diameter), and then inoculated in petri dishes (10 cm diameter) containing 35 mL of potato dextrose agar (PDA) medium as described previously (Gao et al., 2022; Guo et al., 2023). Prior to inoculation, a sterilized cellophane membrane was placed on the surface of the PDA medium for efficient collection of the mycelium samples. The two strains were cultivated in a growth incubator at 25°C and in permanent darkness. After 6 days of growth at 25°C, the fungal cultures were divided into control and treatment groups, in which the latter were subjected to heat exposure at 37°C for 12 h.

Fungal colony diameter was measured every two days once the inoculant was germinated. The inhibition rate was calculated as follows:

where G1 and G2 denote the growth rate of the fungal colony after and before heat treatment. Growth rate was calculated according to the diameter change during the period of growth before and after heat treatment (mm/d).

Reactive oxygen species (ROS) production was measured according to the previously described method (Ren et al., 2017; Guo et al., 2023) with modifications. Briefly, for 2′,7′-dichlorodihydrofluorescein diacetate (DCFH-DA) staining, sterile glass coverslips were obliquely inserted into the petri dishes after the inoculation of fungal blocks, allowing the fungal mycelium to grow later. When hyphae grew on the coverslips, the coverslips with hyphae were incubated in 10 μmol/L of DCHF-DA (Solarbio Life Sciences, Beijing, China) phosphate-buffered saline (PBS) solution for 25 min under dark conditions for ROS visualization. After staining, ROS production was visualized using a fluorescence microscope (Olympus, IX71, Tokyo, Japan). For nitroblue tetrazolium (NBT, Beyotime Biotechnology, Beijing, China) staining, mycelial pellets were incubated in 0.5 μg/mL NBT aqueous solution for 2 h, and the intensity was monitored under a light microscope. The hydrogen peroxide (H₂O₂) content was measured using assay kit following the manufacturer’s instruction (Solarbio Life Sciences, Beijing, China). For measurement of MDA content, 0.1 g of fresh mycelium samples were used for MDA extraction following the manufacture’s instruction (Solarbio Life Sciences, Beijing, China).

Total RNA was extracted from mycelia of L. edodes (five replicates) using TRIzol® Reagent according to the manufacturer’s instruction (Invitrogen, Carlsbad, CA, USA). The RNA quantity determination, libraries construction, and Illumina (Illumina, San Diego, CA, USA) sequencing were performed according to the protocols described by Guo et al. (2023). The raw paired end reads were trimmed, and quality controlled by SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters. Clean reads were separately aligned to reference genome with orientation mode using HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml) software (Kim et al., 2015). The reference genome sequence data reported in this paper have been deposited in the Genome Warehouse (Chen et al., 2021) in National Genomics Data Center, Beijing Institute of Genomics, Chinese Academy of Sciences / China National Center for Bioinformation, under accession number GWHBGWX00000000 that is publicly accessible at https://ngdc.cncb.ac.cn/gwh/Assembly/23974/show. The mapped reads of each sample were assembled by StringTie in a reference-based approach (Pertea et al., 2015).

Gene expression levels (fragments per kilobases per million reads) were calculated using the RESM software (http://deweylab.biostat.wisc.edu/rsem/) (Li and Dewey, 2011). Differentially expressed genes (DEGs) were identified using the DESeq2 (Love et al., 2014) with the following general criteria: |log2FC| > 1 and adjusted p value (padjust) ≤ 0.05. The padjust was calculated using the BH (FDR correction with Benjamini/Hochberg) methods (Benjamini and Hochberg, 1995). Transcripts were annotated using the functional database including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), NCBI non-redundant protein sequences (NR), Swiss-Prot, Protein family (Pfam), and Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups (EggNOG) (Zhang et al., 2024).

Fifty mg of mycelium was used for metabolite measurements using an ultra-high performance liquid chromatography (UHPLC) system (Thermo Electron Corporation, San Jose, CA, USA), coupled with Thermo UHPLC-Q Exactive Mass Spectrometer (MS/MS) (Thermo Electron Corporation, San Jose, CA, USA). Mycelium sample preparation, metabolite extraction, chromatographical measurement, and raw chromatographic data analysis were performed according to the previously described protocol (Guo et al., 2023). Metabolic features detected at least 80% in any set of samples were retained. After filtering, minimum metabolite values were imputed for specific samples in which the metabolite levels fell below the lower limit of quantitation and each metabolic features were normalized by sum. The internal standard (2-Chloro-L-Phenylalanine, HPLC hyper grade, Merck, Darmstadt, Germany, 0.02 mg/mL) was used for quantitative determination of metabolites. Metabolic features which the relative standard deviation (RSD) of QC > 30% were discarded. Compounds were identified by comparison of accurate mass, MS/MS fragments spectra and isotope ratio difference against biochemical databases Human metabolome database (HMDB) (http://www.hmdb.ca/) and Metlin database (https://metlin.scripps.edu/). The mass tolerance between the measured m/z values and the exact mass of the components of interest was ±10 ppm. Mass features without MS/MS spectra were tentatively annotated on MS1 level using 5.0 mDa tolerance. The metabolite abundances were quantified according to their peak areas.

Differentially accumulated metabolites (DAMs) (VIP in OPLSDA model ≥1, false discovery rate (FDR) ≤ 0.05, Fold change ≥ 1) between groups were mapped into biochemical pathways through metabolic enrichment and pathway analysis based on database search (KEGG, http://www.genome.jp/kegg/). To identify statistically significantly enriched pathway, scipy.stats (Python packages) (https://docs.scipy.org/doc/scipy/) was used with Fisher’s exact test.

Principal component analysis (PCA) and orthogonal partial least square regression discriminant analysis (OPLS-DA) were performed to examine overall variations of the metabolomics and transcriptomic data. PCA was conducted using the “FactoMineR” (Lê et al., 2008) and “factoextra” (Kassambara and Mundt, 2017) package in R (version 4.4.1). OPLS-DA was performed using the SIMCA 14.1 (Umetrics, Umeå, Sweden). The robustness of OPLS-DA models was assessed by seven-fold cross-validation. The reliability and accuracy of the predictive models was assessed by analysis of variance testing of cross-validated predictive residuals (CV-ANOVA), R2 and Q2metrics. The p value was estimated with paired student’s t-test on single dimensional statistical analysis. To test the overall associations between DEGs and DAMs, we performed Procrustes analysis using the R package “vegan” (version 2.7–0) (Dixon, 2003). Data were logarithmically (log10) transformed, centered, and pareto scaled prior to multivariate analyses. For integration and visualization of transcriptome × metabolome associations, the sparse partial least squares (sPLS) regression method (Lê Cao et al., 2008) was implemented using the R package “mixOmics” (Rohart et al., 2017).

We first examined the heat stress (HS) responses of the two sexually compatible monokaryons on the morphological and physiological levels (Figure 1). Results demonstrated that HS markedly inhibited the growth rate of mycelium for both monokaryons. When HS was removed, the regenerated mycelium appeared weak compared to the old mycelium formed prior to HS exposure (Figures 1A,B). Based on the average growth rate of the two monokaryons before and after HS, we found that the growth inhibition rate of SP3 strains was significantly higher than that of SP30 strains (Figure 1C). This suggested that the growth of SP3 might be more inhibited than SP30 subjected to HS. DCFH-DA and NBT staining showed that ROS contents were comparable forboth strains at room temperature, whereas were significantly increased when challenged with HS (Figures 1D,F). We also measured the H₂O₂ content by enzymatic assay, which showed consistent results with the ROS staining assay. The MDA production, a marker of the cell membrane lipid peroxidation, was increased significantly under HS (Figure 1E).

Non-targeted metabolomics was applied to examine the metabolic variations of mycelium of the two monokaryons. In total, 1,122 metabolites were detected (Supplementary Table 1). Based on the HMDB superclass, 33.59% of those compounds were lipids and lipid-like molecules, and 23.33% of those were organic acids and derivatives (Supplementary Figure 1A). KEGG functional annotation analysis demonstrated that most of those compounds were related to the category of metabolism, followed by the environmental information processing and genetic information processing (Supplementary Figure 1B).

The upset plot showed that 285 and 382 differentially accumulated metabolites (DAMs) were found for SP3 and SP30 subjected to HS, respectively. The number of the up-regulated metabolites in SP30 was 1.8 times higher than that in SP3. There were 58 and 141 DAMs were specifically up-regulated for SP3 and SP30 under HS, respectively (Figure 2A). PCA revealed that the metabolic profiles of the two strains exhibited distinct differences and were also distinguishable from those exposed to HS. Interestingly, the first component (PC1), which explained 33.0% of the total variation, could explain the differences resulting from the two strains, whereas the second component (PC2), which explained 17.1% of the total variation, could explain the variation derived from the HS treatment (Figure 2B). The loading plot depicted the top 20 compounds contributed to the PCA model (Figure 2C). Figure 2D demonstrated the top 20 metabolites contributed to PC1 and PC2, respectively. With the heatmap analysis we could observe that those metabolites could classify into 3 clusters (a, b, c). The metabolites in cluster “a” were clearly related to the HS condition, whereas the cluster “b” and “c” were associated with monokaryon SP30 and SP3, respectively (Figure 2E).

To identify genes and corresponding functional pathways associated with HS for the two monokaryons, we performed a comparative transcriptomics analysis. A total of 129.13 Gb clean data was obtained with a high score of quality (Q30 rate  >  93.58%) (Supplementary Table 2). The average mapping rate was 89.96%, and more than 84% of genes were coding sequencing for all samples (Supplementary Tables 3, 4). The functional annotations could be found in Supplementary Table 5.

In consistent with the metabolomic data, the transcriptomic profiles of the two monokaryons showed distinct patterns under normal and HS conditions. The first component (PC1), which explained 65.4% of the total variation, could explain the differences resulting from the two strains, while the second component (PC1), which explained 13.44% of the total variation, could explain the variation resulting from the HS treatment (Figure 3A). There were 1,484 and 1899 DEGs identified from SP3 and SP30 in response to HS, respectively. The number of both the up and down-regulated genes in SP30 was higher than that in SP3 in response to HS, particularly for the up-regulated genes in SP30, which were 1.6 times higher than those in SP3 (Figures 3B,C). Using the KEGG enrichment analysis, we observed that the DEGs in SP3 subjected to HS were significantly enriched in environmental information processing (MAPK signalling pathways), cell cycle, and metabolism (biosynthesis of nucleotide sugar, and amino sugar and nucleotide sugar metabolism), whereas those in SP30 were glyoxylate and dicarboxylate metabolism, and genetic information processing (protein processing in endoplasmic reticulum) (Figures 3D,E).

With the volcano analysis, we found that 418 metabolites were differentially accumulated between monokaryon SP3 and SP30 (Figure 4A). Those DAMs were mainly enriched in the pathways of the “glycerophospholipid metabolism,” “alpha-linolenic acid metabolism,” “nucleotide metabolism,” etc. (Figure 4B). There were 124 metabolites that were up-regulated, and 181 metabolites that were down-regulated when SP3 was treated with heat (Figure 4C). These DAMs in SP3 were mainly enriched in the KEGG pathways of “glycerophospholipid metabolism,” “alpha-linolenic acid metabolism,” “biosynthesis of cofactors,” etc. (Figure 4D). As for SP30, 222 metabolites were up-regulated, and 159 metabolites were down-regulated under HS (Figure 4E). These DAMs in SP30 were largely enriched in the KEGG pathways of the “alpha-linolenic acid metabolism,” “linolenic acid metabolism,” “glycerophospholipid metabolism” as shown in Figure 4F. Interestingly, for SP3 most of the top 20 altered metabolic pathways tended to be down-regulated, whereas in SP30 those presented an opposite trend. Of the top 20 altered metabolic pathways, 11 pathways were the same, while the trends were different. It was worth to noting that the top changed pathway “glycerophospholipid metabolism” was significantly (p  <  0.001) down-regulated in SP3, while was significantly up-regulated (p  <  0.01) in SP30, highlighting the strain-specific differences of the metabolic profile in response to HS.

Glycerophospholipid metabolism has been shown to be involved in heat stress responses in both plants (Sun et al., 2020; Feng et al., 2023; Zhou et al., 2024) and edible fungi (Yan et al., 2020; Jiang et al., 2022; Guo et al., 2023). The metabolomic data showed that 16 and 13 DAMs were detected in SP3 and SP30 under HS in the glycerophospholipid metabolism pathway, respectively (Supplementary Tables 6, 7), and the transcriptomics showed that only 3 and 10 DEGs were found in the pathway in SP3 and SP30 under HS, respectively (Supplementary Table 8). According to the KEGG pathway, we mapped those genes involved in the glycerophospholipid metabolism (Figures 5A,B). It showed that the glycerophospholipid metabolism was modulated differently in response to HS in the two monokaryons. Upon exposure to HS, only one gene was involved in the conversion of phosphatidyl -serine to phosphatidylethanolamine in SP3, whereas four genes were involved n SP30. The phosphatidylethanolamine was then catalyzed by different enzymes, which synthesized different downstream metabolites. To discover more DEGs which are putatively associated with regulation of glycerophospholipid metabolism, we conducted an association analysis using the DAMs detected in glycerophospholipid metabolism pathway and all DEGs. The arrow plot demonstrated a high agreement between the DAMs and DEGs datasets and the discriminative ability of each dataset (Supplementary Figure 2). Using the sPLS model, we were able to detect two highly associated communities (correlation coefficients ≥ 0.9 ) in the two monokaryons, respectively (Figure 5). The communities of SP3 consisted of 11 genes and 3 metabolites, with 2 metabolites were correlated with 7 genes in subcommunity I, and 1 metabolite was correlated with 4 genes in subcommunity II. As for SP30, 2 metabolites were correlated with 4 genes in subcommunity I, and 2 metabolites were correlated with 3 genes in subcommunity II (Figures 5C–F). For SP3, it was worth noticing that most of the highly associated genes were annotated as being involved in MAPK signaling pathways, which were completely different from what was observed in SP30 (Figures 5E,F).

To unearth the gene regulatory network of the upregulated metabolites elicited by HS, an integrative analysis was performed using the upregulated metabolites with the DEGs of the two monokaryons. First, we examined the agreement between the upregulated metabolites and DEGs data sets using the sPLS regression model. The resulting arrow plots demonstrated a good association between the two datasets for both monokaryons (Figures 6A,B). All the identified highly associated DAMs and DEGs are presented in the Supplementary Figure 9. For strain SP3, the integration analysis revealed two distinct regulatory subnetworks in response to HS. The large subnetwork (subnetwork I) comprised 33 metabolites and 8 putative regulatory genes, and the subnetwork II comprised 4 metabolites and 3 putative regulatory genes. Interestingly, the metabolites in subnetwork I were all negatively correlated with genes, whereas in subnetwork II they were all positively correlated. The metabolites in subnetwork II were negatively correlated with genes in subnetwork I (Figure 6C). The top 5 hub genes included pb011955, pb003418, gb002771, pb007688, gb006174, and 3 of those genes were annotated to heat shock protein (Table 1). For SP30, 2 separated subnetworks were found, of which the smaller subnetworks were partitioned to 2 subcommunities (subnetwork II and subnetwork III). The subnetwork I included negatively correlated metabolites and DEGs, whereas the subnetworks II and III consisted of positively correlated metabolites and DEGs (Figure 6D). The top 5 hub genes were pb006533, pb003154, gb005258, gb007318, pb008404 in SP30, which were completely different from what discovered in SP3 (Table 1).