Transcriptome data reveal conserved patterns of fruiting body development and response to heat stress in the mushroom-forming fungus Flammulina filiformis

Mushroom-forming fungi are complex multicellular organisms that form the basis of a large industry, yet, our understanding of the mechanisms of mushroom development and its responses to various stresses remains limited. The winter mushroom (Flammulina filiformis) is cultivated at a large commercial scale in East Asia and is a species with a preference for low temperatures. This study investigated fruiting body development in F. filiformis by comparing transcriptomes of 4 developmental stages, and compared the developmental genes to a 200-genome dataset to identify conserved genes involved in fruiting body development, and examined the response of heat sensitive and -resistant strains to heat stress. Our data revealed widely conserved genes involved in primordium development of F. filiformis, many of which originated before the emergence of the Agaricomycetes, indicating co-option for complex multicellularity during evolution. We also revealed several notable fruiting-specific genes, including the genes with conserved stipe-specific expression patterns and the others which related to sexual development, water absorption, basidium formation and sporulation, among others. Comparative analysis revealed that heat stress induced more genes in the heat resistant strain (M1) than in the heat sensitive one (XR). Of particular importance are the hsp70, hsp90 and fes1 genes, which may facilitate the adjustment to heat stress in the early stages of fruiting body development. These data highlighted novel genes involved in complex multicellular development in fungi and aid further studies on gene function and efforts to improve the productivity and heat tolerance in mushroom-forming fungi.

Introduction development of F. filiformis or its closely related F. velutipes, our knowledge about the fruiting body development and heat stress response of this mushroom is still incomplete.
In this study, we aimed to (i) uncover key fruiting body genes in various developmental stages; (ii) investigate whether conserved developmental patterns exist in F. filiformis and other Agaricomycetes; and (iii) understand the responses to heat stress of heat sensitive (XR) and resistant (M1) strains and identify the key heat stress response genes in F. filiformis. We sampled RNA in different developmental stages of the M1 strain, and from the primordium stage grown at 10˚C and 18˚C of both the M1 and XR strains, for studying fruiting body development and its response to heat stress, respectively. Our results identified conserved gene expression patterns of fruiting body development in the Agaricomycetes and revealed that the heat tolerant strain M1 differentially expressed more genes in response to heat stress than the heat sensitive strain XR.

Strains and culture conditions
The heat tolerant strain M1 (CGMCC5.2219) was domesticated from a wild strain collected in subtropical areas in China in summer (Fujian province). The heat sensitive strain XR (CGMCC5.2218) was isolated from a mushroom market, this strain was imported from Japan [33]. Both of them are deposit in the Chinese General Microbiological Culture Collection Center (CGMCC). They were grown on enriched Potato Dextrose Agar (PDA) medium (0.05% KH 2 PO 4 , 0.05% MgSO 4 , 2% glucose, 0.2% yeast extract, 0.2% peptone, 1.8% agar) in 90 mm Petri dish for 10 days at 23˚C. Then, the mycelium was inoculated in liquid cultures in 500ml Erlenmeyer flask containing 200ml enriched Potato Dextrose Broth (PDB) medium (0.05% KH 2 PO 4 , 0.05% MgSO 4 , 2% glucose, 0.2% yeast extract, 0.2% peptone), shaken at 150 r.p.m for 10 days at 23˚C. Afterwards, the mycelium was inoculated to a growth medium consisting of 90% cottonseed hull, 10% wheat bran, and 65% water in 1100 ml disposable bottles after sterilization. Inoculated bottles were incubated at 23˚C under dark conditions for 30 days, and then the mycelium scratched to emulates disturbance and transferred to 23˚C with 95% humidity for fruiting.

Sample collection for RNA-seq
We selected different developmental stages of strain M1 grown at 10˚C (except for the mycelium, which was grown at 23˚C). We collected the vegetative mycelium (VM), the primordium (P10-M1), the young fruiting body cap (YFBC), the fruiting body cap (FBC), the young fruiting body stipe (YFBS), and the fruiting body stipe (FBS) (Fig 1A). For the heat stress response study, we collected the primordia, young fruiting body, and fruiting body of M1 strain grown at 18˚C (P18-M1), as well as the primordia, young fruiting body and fruiting body of XR strain grown at 10˚C and 18˚C respectively (Fig 1B).

Total RNA preparation and transcriptome sequencing
Samples from three biological replicates were flash-frozen in liquid nitrogen and stored at -80 C. Total RNA of each sample was extracted using the RNAprep Pure toolkit, following the manufacturer's protocol (TIANGEN, Beijing, China). Sequencing libraries were generated using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's recommendations and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated. All raw-sequence reads data were deposited in NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/ Traces/sra) with accession number of PRJNA557510.

Read mapping to the reference genome, FPKM and gene annotation
Clean data were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads from raw data through Trimmomatic (v.0.33) [41]. Then, we used HISAT2 (v.2.10) [42] to map the clean reads to the reference genome which assembled to chromosome level [20], and employed StringTie (v.1.3.4) [43] to calculate each gene's FPKM value. All of the genes were annotated using local BLASTX programs against the Nr, SwissPort, GO and PFAM databases.

Differential gene expression analysis
We performed differential gene expression analyses on each adjacent developmental stage of the M1 strain ( Fig 1A) as well as the primordium stage grown at 10˚C vs 18˚C of the M1 and XR strains (Fig 1B). Analyses were based on three biological replicates per condition and were performed using the DESeq package (1.18.0) [46]. The input data of young fruiting body grown at 10˚C of M1 (YFB10-M1) were used as the pooled reads of young fruiting body cap (YFBC) and young fruiting body stipe (YFBS) of M1, and the input data of fruiting body grown at 10˚C of M1 (FB10-M1) were used as the pooled reads of fruiting body cap (FBC) and fruiting stipe (FBS) of M1. Genes with log2 (fold change) � 1 and Padj � 0.05 were considered as differentially expressed gene (DEG).

Gene network construction and visualization
Co-expression networks were constructed using the WGCNA package in R [47]. Genes with averaged FPKM from three replicates higher than 1 in at least one sample were input to WGCNA unsigned co-expression network analysis (S1 Table). The modules were obtained using the step-by-step network construction function on block-wise modules with default settings, except that the power is 24 for fruiting body development analysis, 20 for M1 and 10 for XR in heat stress response analysis. TOM-Type was set to signed, minModuleSize to 30, and mergeCutHeight to 0.25. The networks were visualized using Cytoscape (v.3.5.1) [48].

Comparative genomic analysis
In order to check the conservation level of the developmentally regulated genes? during fruiting body development, analyzed a 201 genome dataset (ranging from unicellular yeasts to filamentous and complex multicellular fungi which also included F. filiformis in this study) and the corresponding phylogenetic tree taken from a previous study [13]. Conservation of genes was assessed based on the phylogeny, by assessing the presence/absence of genes across the panel of species.

qRT-PCR analysis
Reverse transcription of RNA (1ug) in a 20 μL reaction volume was performed using TUREscript 1st Stand cDNA SYNTHESIS Kit, following the manufacturer's protocol (Aidlab, Beijing, China). Reactions were incubated at 42˚C for 60min, and 65˚C for 10min. The amplifications were performed using 5 μL SYBR qPCR Mix, 0.5 μL forward primer, 0.5 μL reverse primer and 1 μL cDNA, and 3 μL ddH2O in a final volume of 10 μL. The cycling parameters were 95˚C for 3 min followed by 30 cycles of 95˚C for 10 s, 58˚C for 30 s and 72˚C for 20 s. The relative gene expression was analyzed calculated by the qPCRsoft3.2. The 18S ribosomal RNA gene was used as the internal reference. The primers of each gene were listed in S2 Table.

Overview of the transcriptome sequence data
We obtained 39.4-63.6 million paired-end reads for 15 sample types in triplicates (45 libraries in total) (S3 Table). Of the quality-filtered reads 58.3-71.2% mapped to the reference genome of F. filiformis (S3 Table). The moderate mapability may be caused by the different strains used in this study compared with the reference genome from strain KACC42780 [20]. Although we expect this to not influence our results, it may cause lower sensitivity in faster evolving or strain specific genes. Differences in read mapability were not found between the M1 and XR strain which made the transcriptome comparable across these two strains (S3 Table). To validate the results of the RNA-seq analysis, 18 genes were randomly selected for quantitative realtime PCR (qRT-PCR). These genes showed expression patterns similar to those in the RNAseq data (S1 Fig), indicating that our transcriptome sequencing provided a good estimate of gene expression patterns in the analyses of fruiting body development and heat stress response of F. filiformis.  Table). This gene expression pattern was also recognized by the WGCNA analysis, a systems biology approach aimed at uncover gene modules which share gene expression patterns at a pre-specified similarity cutoff [47]. We identified six gene modules highly correlated with a single tissue type (Fig 2A; S5 Table). Among them, the primordium module (module no. 6), young fruiting body cap module (module no. 4) and mature fruiting body cap module (module no. 12) contained the highest number of genes (Fig 2A). These results indicate that primordium stage comprises the most significant morphogenetic transition, and that hymenium maturation and sporulation in young and mature fruiting body caps may also harbor complex molecular mechanisms. This gene expression pattern is consistent with those found in other Agaricomycetes [8,22,13,49]. Moreover, the DEGs related to each developmental stage were enriched for GO terms typical for fruiting body formation, see S3

Primordium development includes genes widely conserved in Agaricomycetes
The top 20 up-regulated genes induced in primordium relative to vegetative mycelium were listed in Table 1. Building on a previously published dataset [13], we found these genes were widely conserved in fungi, and re-emphasized that several primordium-upregulated genes have homologs in simple multicellular or unicellular fungi. This indicates that some conserved gene families were recruited for complex multicellularity during evolution [13]. Among them, the Gti1/Pac2 family is conserved in all fungi (Table 1; Fig 3). Previously, this family has been discussed mostly in yeasts and pathogenic fungi, where it plays an important role in fungal growth and development [50]. Recent research revealed that this family is also developmentally regulated in Armillaria ostoyae, Coprinopsis cinerea, Lentinus tigrinus, Rickenella mellea, Schizophyllum commune, and Phanerochaete chrysosporium [13]. Thus, these genes may also play a key role in fruiting body development in Agaricomycetes. In addition, three TFs in the top 20 up-regulated genes (Zinc finger, C2H2 type, Zinc finger, Ring type, and Zn (2)-C6 fungal type) are conserved in Dikarya and Zoopagomycota plus later diverging phyla (Table 1; Fig 3), which reinforces the role of these TFs in complex multicellularity in fungi [13].
We found the gene encoding flammutoxin was conserved in 19 species in Agaricomycetes and 1 species in the Dacrymycetes (Table 1; Fig 3). This protein has been studied in F. filiformis, and may form a pore in the intestinal epithelial cells of fungivorous animals, leading to cell death [51]. However, because this protein is heat-labile, clinical reports about the intestinal dysfunction caused by ingestion of this mushroom are rare. Pore-forming proteins have been studied in Pleurotus species, which exhibit cytotoxicity toward insect cells via pore formation in cell membranes to defend predation [52,53]. Thus, whether flammutoxin serves to protect the fruiting body of F. filiformis predated by mammals or insects needs further investigation. In addition, a ricin-B lectin gene (chr11_AA_01461), showed high expression in fruiting bodies. It is homologous to the Macrolepiota procera mpl, Mpl protein shows toxicity towards the nematode Caenorhabditis elegans [54]. A re-analysis of its homologs in C. cinerea, L. tigrinus, R. mellea, S. commune, and P. chrysosporium revealed that it has fruiting body specific expression pattern in these species as well. This is in concordance with its expression pattern in M. procera, indicating a conserved role of this gene in protecting fruiting bodies against predators and parasites in Agaricomycetes [54]. Compared with the conserved flammutoxin and ricin-B lectin gene discussed here, we found the previously mentioned widely conserved ribotoxins gene is not exists in F. filiformis genome [55].
Interestingly, the second most up-regulated gene in primordium was restricted to the family Physalacriaceae, which indicated a role in primordium formation in this family (Table 1;  Fig 3). Unfortunately, our knowledge about this gene is limited. However, recent study reported the successful overexpression and RNA interference of the transcription factor pdd1 in F. filiformis, which provides genetic tools to study this gene in the future [56]. Based on our analysis, almost half of the top 20 up-regulated genes in primordium stage without annotations (Table 1), they are widely conserved in simple multicellular and complex multicellular fungi, indicating that systematic studies are needed on mushroom development to help to understand multicellularity.

Notable genes during fruiting body development
In the primordium stage, we found that genes encoding hydrophobins were (homologs of fv-hyd1, fv-hyd7 and S. commune sc3) significantly up-regulated relative to vegetative mycelium (Fig 3; S2A Fig; S4 Table). It is well known that hydrophobins and cerato-platanins assemble at the hyphae surface to promote their aggregation in humid environments [57][58][59]. Hydrophobins may also hinder water absorption through the membrane [60]. We detected the aquaporin gene specifically induced in primordia and stipe tissues (Fig 2B;  S2A Fig; S4 Table). Aquaporins are integral membrane proteins responsible for water and solute transport, and also involved in mycorrhizal formation and plant-fungal interactions during symbiosis establishment [26,61]. Recent study revealed that aquaporins were also developmentally regulated in L. bicolor fruiting bodies [26]. However, aquaporins have been discussed mostly in the context of mycorrhizal fungi [62][63][64][65]. Based on our results, it seems likely that aquaporindependent water transport is a key process during mushroom development in saprotrophic fungi too, possibly in water transport along the stipe to facilitate water supply of the developing cap and gills [64]. A re-analysis of the data published by Sipos et al. (2017) [22] and Krizsán et al. (2019) [13] indicated that aquaporins were developmentally regulated in all six species studied by these authors, indicating that the role of aquaporins in development is not restricted to Flammulina, but may be widely conserved in fruiting body development. Aquaporins were differentially expressed in mycorrhizal species [2,26,64], which provides additional support to the hypothesis that fruiting bodies and ectomycorrhizae have many shared gene expression patterns, possibly pointing to common developmental origins [2,64,66].
Interestingly, we found that the Flammulina homolog of the stipe elongation gene eln3 of C. cinerea possessed similar expression pattern to those of aquaporin genes during fruiting body development (Figs 2B and 3; S2A Fig; S4 Table). The mutant strain of this gene in C. cinerea produced aberrant fruiting bodies, in which the stipe hardly elongated during development [67]. The eln3 homolog of Volvariella volvacea was also reported to be differentially expressed during fruiting body development [68]. A re-analysis indicated that eln3 in C. cinerea and its homologs in L. tigrinus, A. ostoyae and R. mellea were developmentally regulated in RNA-Seq on data from previous studies [22,13]. Homologs in S. commune and P. chrysosporium were not developmentally expressed, which might be explained by the lack of a stipe in these species. The broad conservation and expression patterns of eln3 suggests that the molecular mechanisms of stipe elongation may be shared in Agaricomycetes, despite the independent origins of pileate-stipitate fruiting bodies in the class [28]. These results further highlight this gene as an interesting target in future strain improvement programs.
We found that one of the hub genes in the primordium module (module no. 6, WGCNA co-expression analysis) was a homolog of Aspergillus nidulans ppoA (Figs 2C and 3). This gene participates in oxylipin synthesis, which modulates sexual and asexual development in A. nidulans [69]. During sexual development, the PpoA protein initially localized in Hülle cells formed at the stage of cleistothecial primordium formation, and subsequently in immature cleistothecia in A. nidulans [69]. Over-expression of this gene in A. nidulans promotes sexual spore formation [69]. A re-analysis of this gene's homologs in C. cinerea, R. mellea, and A. ostoyae revealed a conserved expression pattern in these species, which indicates that oxylipins may mediate sexual development in the Ascomycota and the Agaricomycetes.
The homolog of Lentinula edodes priA was highly induced in primordium and young fruiting body cap of F. filiformis (S2A Fig; S4 Table). This gene was reported to possess the higher expression level in primordium and young fruiting body of L. edodes and over-expressing the priA gene in L. edodes monokaryotic mycelium remarkably decreased zinc ion accumulation, which indicates this gene may play a role in regulation of the intracellular zinc concentration [70]. Surprisingly, we found its homologs were highly expressed in vegetative mycelium and lower expressed in fruiting bodies in A. ostoyae, C. cinerea, L. tigrinus, R. mellea, S. commune, and P. chrysosporium.
A large number (74) of ribosomal protein encoding genes, and homologs of genes involved in cell differentiation and cell wall formation in S. pombe (agn1 and rae1), Candida albicans (ada2), and S. cerevisiae (ecm3) were hub genes in the young fruiting body cap module of the WGCNA analysis or were up-regulated relative to the primordium stage (Figs 2D and 3; S2A and S2B Fig; S4 Table). These results might reflect intense growth and protein synthesis in young fruiting bodies. Although the hymenium was immature in this stage, four meiosis regulation genes homologous to S. pombe mei2, Cryptococcus neoformans dmc1, and S. cerevisiae hfm1, hop1, and zip4 were up-regulated (Fig 3; S2A Fig; S4 Table). Among them, dmc1, hop1 and hfm1 homologs were also up-regulated in mature cap (Fig 3; S2A Fig; S4 Table).
Compared to young fruiting body cap, more meiosis genes (homologs of S. pombe psm1 gene, C. neoformans dmc1 gene, and S. cerevisiae, msh4, msh5, meu13, hop1, hfm1, and smg7) were induced in the fruiting body cap (Fig 3; S2A Fig; S4 Table). Among the genes induced in young fruiting body cap and fruiting body cap, the genes homologous to C. neoformans pum1 may be noteworthy (Fig 3; S2A Fig; S4 Table). Pum1 is an RNA binding protein, and possesses an important role in post-transcriptional regulation in basidium development and sporulation in C. neoformans [71][72][73]. Previous studies revealed that the knockout of this gene in C. neoformans resulted in a severe defect in basidium formation [71][72][73]. In this study, we detected five genes homologous to C. neoformans pum1, two of them were developmentally regulated in young fruiting body cap and fruiting body cap, which indicates they may participate in basidium formation and sporulation in F. filiformis (Fig 3; S2A Fig). We detected another sporulation-related gene, homologous to A. nidulans samB (Fig 3; S2A Fig). Knock out of this gene in A. nidulans hindered ascospore formation [74]. A re-analysis of the homologs of pum1 and Samb in A. ostoyae, C. cinerea, L. tigrinus, R. mellea, S. commune, and P. chrysosporium revealed they possess conserved expression patterns in these species. Due to these two genes were widely conserved (Fig 3), we therefore speculate that some molecular mechanisms of spore formation may be conserved in fungi.

CAZymes and F-box genes
Certain CAZymes were shown or assumed to participate in cell wall remodeling during fungal tissue differentiation [13,49,[75][76][77]. We annotated 407 CAZymes genes in F. filiformis, 137 of them were differentially expressed (Fig 2E; S7 Table), which is consistent with previous studies in other mushroom-forming fungi [13,20,49]. Among these genes, Glycoside hydrolases (GH) and Glycosyltransferases (GT) were most abundant, with 57 and 30 genes, respectively. Although the targets of these families in fruiting bodies are currently unknown, the stage specific expression of these genes during fruiting body development reinforces the view that cell wall remodeling is a widespread and well-organized process in fruiting body development in Agaricomycetes.
F-box proteins play a key role in protein ubiquitination and modification, and are involved in many important biological processes not only in plants, but all eukaryotic [78,79]. They were recently reported in relation to fruiting body development [13,49]. In this study, 210 Fbox encoding genes were annotated in F. filiformis, of which 80 were developmentally expressed and showed stage-specific expression patterns (Fig 2E; S8 Table). Similar expression patterns were also recognized in A. ostoyae, C. cinerea, R. mellea, L. tigrinus, and S. commune, which suggests that F-box genes may be crucial during fruiting body development in Agaricomycetes [13].

Strain M1 expressed a large gene pool in response to heat stress relative to XR
Based on cultivation tests, we found that the growth of M1 strain showed no difference in 10˚C and 18˚C, while, the growth in XR strain was obviously retarded at 18˚C. On the molecular level, we found more differentially expressed genes in M1 (882 DEGs) than in XR (556 DEGs) (Fig 1C; S4 Fig). Based on the Venn diagram on Fig 1, we found two genes with elevated expression level in all developmental stages of strain M1 at 18˚C: a DNA damage repair gene homologous to S. cerevisiae rad18 and a epoxide hydrolase gene homologous to A. niger. For strain XR, seven genes had an elevated expression level at all developmental stages under 18˚C (Fig 1C): two hsp20 genes, one WD repeat-containing gene, and four genes without annotation. Consistent with these functions, a GO enrichment analysis revealed that up-regu-  Table). These results indicate the different heat stress response strategies were employed in these two strains.
Specifically, differential expression analyses revealed that the M1 strain had more heat shock protein genes up-regulated than XR (S4 Fig). Among them, homologs of S. cerevisiae hsp70, S. pombe hsp90 and another gene homologous to Ustilago maydis fes1 induced in M1 at 18˚C may be noteworthy. The Hsp70 protein could protect nascent polypeptides and refold the damaged proteins under heat stress conditions [80]. If protein folding fails with Hsp70, Fes1 could interact with misfolded proteins and lead to their destruction by the ubiquitinproteasome machinery [81]. Compared with Hsp70, Hsp90 functions primarily in the final maturation of proteins. Therefore, these genes may act as an "assembly line" [80] of protein maturation under heat stress during primordium development of strain M1. The heat stress induction of hsp70 and hsp90 was also reported in Lentinula edodes and Ganoderma lucidum [82,83]. Homologs of these two genes were not differentially expressed in XR strain. Instead, the genes homologous to S. pombe hsp20 were up-regulated in XR strain in all developmental stages at 18˚C, they play different roles than hsp70 and hsp90, which are probably required to prevent misfolded protein aggregation and their degradation under heat stress [80].
In addition, the DNA damage repair gene homologous to S. cerevisiae rad18 was up-regulated in all developmental stages in M1. Yeast strains lacking Rad18 proteins may be highly sensitive to a wide variety of DNA damaging agents such as UVC-light, ROS stress and γ-radiation [84][85][86]. However, this gene was only up-regulated in the young fruiting body stage of XR grown at 18˚C. This example, combined with the fewer DEGs in the XR strain may indicate the loss of an ancestral heat stress response mechanisms in the commercial XR strain. This may result from the lack of temperature fluctuations and stress in general in factory settings.

Conclusions
This study broadened our knowledge of fruiting body development and heat stress response of mushroom-forming fungi based on comparisons of transcriptomic data in F. filiformis. We detected a series of genes (e.g. aquaporins, eln3-homologs, hydrophobins, conserved transcription factors, oxylipin biosynthesis genes) that show conserved, dynamic expression during fruiting body development, and also uncovered signal for defense against high temperature in the heat tolerance strain (M1) (e.g. hsp70, hsp90 and fes1 homologs). These, or other differentially expressed genes, might be good candidates for in-depth experimental follow-up analyses (e.g. gene knockout) to understand their specific roles and answer important or interesting questions that remained open. Analyzing the function of conserved genes in model and nonmodel species will be necessary to broaden our knowledge on fruiting body development in the Agaricomycetes.