Jia-Ning Wan1†Yan Li1†Ting Guo1Guang-Yan Ji2Shun-Zhen Luo2Kai-Ping Ji2*Yang Cao2
Qi Tan1
Da-Peng Bao1
Rui-Heng Yang1*- 1Key Laboratory of Agricultural Genetics and Breeding of Shanghai, Key Laboratory of Edible Fungal Resources and Utilization (South), National Engineering Research Center of Edible Fungi, Institute of Edible Fungi, Shanghai Academy of Agricultural Sciences, Shanghai, China
- 2Hongzhen Agricultural Science and Technology Co. Ltd., Jinghong, China
Phlebopus portentosus (Berk. and Broome) Boedijin, a widely consumed mushroom in China and Thailand, is the first species in the order Boletaceae to have been industrially cultivated on a large scale. However, to date, the lignocellulose degradation system and molecular basis of fruiting body development in P. portentosus have remained cryptic. In the present study, genome and transcriptome sequencing of P. portentosus was performed during the mycelium (S), primordium (P), and fruiting body (F) stages. A genome of 32.74 Mb with a 48.92% GC content across 62 scaffolds was obtained. A total of 9,464 putative genes were predicted from the genome, of which the number of genes related to plant cell wall-degrading enzymes was much lower than that of some saprophytic mushrooms with specific ectomycorrhizal niches. Principal component analysis of RNA-Seq data revealed that the gene expression profiles at all three stages were different. The low expression of plant cell wall-degrading genes also confirmed the limited ability to degrade lignocellulose. The expression profiles also revealed that some conserved and specific pathways were enriched in the different developmental stages of P. portentosus. Starch and sucrose metabolic pathways were enriched in the mycelium stage, while DNA replication, the proteasome and MAPK signaling pathways may be associated with maturation. These results provide a new perspective for understanding the key pathways and hub genes involved in P. portentosus development.
Introduction
Phlebopus portentosus (Berk. and Broome) Boedijin is an ectomycorrhizal edible mushroom with a widespread distribution in tropical parts of China, especially Yunnan, Guangxi, and Hainan provinces (Zhang et al., 2017). Because it is rich in nutrients, e.g., polysaccharides, amino acids, mineral elements, and pyrrole alkaloids, this fungus is popular and widely used in these regions (Sanmee et al., 2010; Sun et al., 2018; Kumla et al., 2021). P. portentosus has been successfully artificially cultivated and produces sporocarps in artificial substances in vitro in greenhouses or factories in China and Thailand, making it the first species in the order Boletaceae to have been industrially cultivated on a large scale (Ji et al., 2011; Kumla et al., 2015; Zhang et al., 2017). Especially in China, the industrialized cultivation of P. portentosus is highly similar to the Agaricus mushroom industry, with a daily production of 6 tons (Zhang et al., 2017). Nutrition and management techniques play important roles in maintaining the yields and quality of mushrooms (Tseng and Luong, 1984). Furthermore, understanding the lignocellulose degradation system and fruiting body development of this fungus is crucial to enhance yields and improving quality (Isikhuemhen and Mikiashvilli, 2009). However, the molecular basis of these biological processes is not fully understood.
Although P. portentosus can be cultivated by providing carbon or nitrogen compounds in artificial media, it is still characterized as an ectomycorrhizal (ECM) fungus that requires the use of coculture and isotopic methods (Pham et al., 2012; Kumla et al., 2016). Because most ectomycorrhizal fungi have difficulty producing macroscopic sporocarps in axenic culture on a large scale (Hall et al., 2003), some studies have revealed that two ecological niches may be required for P. portentosus growth, both independently as a saprophyte and in association with plants as an ectomycorrhizal symbiont, offering an excellent model to study the mechanisms by which fungi obtain nutrition from in vitro environments (Cao et al., 2015; Zhang et al., 2017). Plant cell wall-degrading enzymes (PCWDEs) secreted by fungi, including cellulases, hemicellulases, pectinases, and ligninases, provide fungi with the means to acquire energy and nutrients from plant cell walls (Hammel, 1997). However, genomic and transcriptomic comparisons between ectomycorrhizal and saprophytic fungi have revealed a reduced repertoire of genes encoding PCWDEs in ectomycorrhizal fungi, resulting in a limited capacity to decompose lignocellulose (Kohler et al., 2015; Miyauchi et al., 2020). A genome analysis of P. portentosus also revealed the presence of few or no copies of carbohydrate metabolism enzymes (CAZymes) that act on cellulose, xylan, pectin, and lignin (Cao et al., 2015). It is surprising that so few genes for lignocellulose degradation support P. portentosus growth on sawdust substances. However, our understanding of these enzymatic repertoires of P. portentosus is far from complete.
The hyphal differentiation, fruiting body formation and development of edible mushrooms are complex processes. The expression of genes at different stages of fungal development (hyphae, primordia, and fruiting body) is associated with yield and quality. Thus, an understanding of the molecular mechanisms regulating fruiting body development is necessary to generate improved strains and varieties of mushrooms by genetic modification or breeding. With the advent of next-generation sequencing techniques, genome and transcriptome sequencing of saprophytic and ectomycorrhizal mushrooms has been widely used to assess gene expression in different growth stages, including for Agrocybe aegerita (Wang et al., 2013), Lentinula edodes (Song et al., 2018; Wang et al., 2018), Flammulina velutipes (Liu et al., 2020), Pleurotus eryngii (Xie et al., 2018), Morchella (Hao et al., 2019; Wei et al., 2019), and Tricholoma matsutake (Tang et al., 2020). Stage-specific pathways enriched and correlated with the growth and development of fruiting bodies have been identified, with functions related to fungal cell wall remodeling, targeted protein degradation, signal transduction, adhesion, and small secreted proteins. Although P. portentosus has been successfully artificially cultivated, gene expression in this fungus grown in vitro during different developmental stages is poorly understood.
To better understand the molecular mechanism of lignocellulose degradation and fruiting body development, whole-genome and transcriptome sequencing of P. portentosus was performed. Three developmental stages were analyzed, including hyphal, primordia, and fruit body stages. The results of genomic and transcriptomic analyses will improve our understanding of the stage-specific expression of genes associated with the functional properties of fruiting body development in P. portentosus. In addition, the transcriptomic profiles of ectomycorrhizal (ECM) fungi cultivated in vitro provide a foundation for further research on the cultivation of other ectomycorrhizal (ECM) fungi.
Materials and Methods
Strain Collection, Artificial Cultivation and Monokaryotic Strain Preparation
The strain PP17026, which was collected from Yunnan Province, China, was provided by Hongzhen Agricultural Science and Technology Co. Ltd. and cultured on M1 medium in the dark at 30°C. Spawn preparation, artificial cultivation and management of P. portentosus was performed according to previously published methods (Yang et al., 2019). Samples were collected in the hyphal (S, stage III), primordia (P, stage VI), and fruiting body (F, stage VII) stages (Yang et al., 2019). The hyphae fully covered on the surface of substrates, the small primordia with about 0.2 cm height and the stipes and pileus (cut into 0.2 cm × 0.2 cm pieces) used as representative samples of S, P, and F stages respectively were collected and flash-frozen in liquid nitrogen and stored at –80°C before RNA extraction.
The monokaryotic protoplasts were isolated from PP17026 using a previously published method (Chang et al., 1985). After 10 days of incubation, mycelia were collected and stored at –80°C until DNA and RNA extraction for genome and transcriptome sequencing. Transcriptome sequencing of monokaryotic strains was performed to correct the gene prediction results from the genome annotation.
DNA Extraction, Library Construction and Genome Sequencing
DNA was extracted from the monokaryotic strain using the CTAB method. A total of 100–200 mg wet weight mycelium was ground in liquid nitrogen using a pestle in a centrifugal tube. Then, 1.5 mL of extraction buffer was added containing 20 mM EDTA, 100 mM Tris-HCl, 1.5 M NaCl, 2% CTAB, and 1% β-mercaptoethanol, and the sample was incubated at 65°C for 30 min. Protein and polysaccharide removal and DNA precipitation were conducted using chloroform-isoamyl alcohol (24:1, v/v) and isopropanol, respectively. The concentration, purity and integrity of DNA were assessed by NanoDrop (Thermo Scientific), Qubit and pulsed field electrophoresis analysis, respectively. A 20-kb sequencing library was built using large DNA segments with an ONT Template prep kit (SQK-LSK109) and an NEB Next FFPE DNA Repair Mix kit. The high-quality library was sequenced on an ONT PromethION platform with a corresponding R9 cell and an ONT sequencing reagent kit (EXP-FLP001.PRO.6). In addition, a small 300-bp sequencing library was built for the Illumina platform to improve the accuracy of the long-fragment library sequencing results.
RNA Extraction, Library Construction and RNA Sequencing
Total RNA were extracted from all samples as previously described using TRIzol reagent (Invitrogen, Burlington, ON, Canada). Three biological replicates were processed. The purity and concentration of RNA was assessed using a NanoPhotometer® spectrophotometer (IMPLEN, CA, United States) and a Qubit® RNA Assay Kit with a Qubit®2.0 Flurometer (Life Technologies, CA, United States). RNA integrity was assessed using an RNA Nano 6000 Assay Kit with the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). Subsequently, RNA was treated with DNase I (Thermo Fisher, Waltham, MA, United States), and the mRNA was purified based on PolyA selection and fragmentation. First strand cDNA synthesis was performed using SuperScript II (Thermo Fisher) followed by second strand cDNA synthesis, end repair, 30-end adenylation, adapter ligation and PCR amplification. Purification was performed using AmPureXP Beads (Beckman Coulter, Brea, CA, United States). The 300 bp libraries were sequenced for paired-end 150-bp-reads on Illumina HiSeq Novaseq 6000 platform (Illumina Inc., United States).
Genome Assembly and Annotations
To retrieve the nucleotide sequences from raw signal data generated from the ONT platform, base calling was performed using Albacore implemented in MinKNOW (Sahoo, 2017). Subsequently, by filtering the low-quality reads and demultiplexing the ONT barcodes and short reads, the filtered subreads were obtained using Canu v1.5 (Sergey et al., 2017). Raw Illumina data trimming was performed according to previously published methods (Yang et al., 2016). The corrected subreads were assembled into contigs using wtdbg v1.2.8 (Ruan and Li, 2020). To improve assembly, contig correction with Illumina read data was performed using Pilon (Walker et al., 2014). The completeness of the genome assembly was evaluated using BUSCO 4.0 with fungi_odb9 (Simão et al., 2015). MicroRNAs, rRNAs, and tRNAs were identified using tRNAscan-SE (Lowe and Eddy, 1997) and Infernal 1.1 (Nawrocki, 2014).
Three methods were used for gene prediction: (i) using a combination of the results from Genscan (Burge and Karlin, 1998), Augustus v2.4 (Mario and Burkhard, 2005), GlimmerHMM v3.0.4 (Majoros et al., 2004), GeneID v1.4 (Blanco et al., 2002), and SNAP (version 2006-07-28) (Ian, 2004); (ii) homologous protein prediction conducted using GeMoMa v1.3.1 (Jens et al., 2019); and (iii) unigenes were assembled and predicted using Hisat2 v2.0.4 (Kim et al., 2019), StringTie v1.2.3 (Pertea et al., 2015) and TransDecoder v2.0 based on the transcriptome sequencing of monokaryotic strains. All three results were integrated using EVM v1.1.1 to obtain an accurate prediction (Haas et al., 2008). Predicted genes were annotated using BLAST (Basic Local Alignment Search Tool) searches against Swiss-Prot (Boeckmann et al., 2005), TrEMBL (Rédei, 2008) and Nr databases. Gene ontology and KEGG metabolic pathway matches were identified using local Blast2GO tools (Conesa et al., 2005) and KAAS (Yuki et al., 2007), respectively. All predicted protein families were analyzed with InterProScan (Quevillon et al., 2005) and Pfam analysis (Finn et al., 2010). Carbohydrate-active enzymes were annotated using dbCAN (Huang et al., 2018). Oxidoreductases were extracted from the proteins predicted from the genome using a combination of IPR domain searches and the JGI cluster pipeline1 according to a previously published method (Floudas et al., 2012). The A- and B-mating-type genes were identified using previously published methods (Chen et al., 2016). The genome has been submitted to the NCBI under the accession numbers of JAHRGP000000000.
Data Processing, Alignment With the Reference Genome, Differentially Expressed Gene Analysis and Annotations for RNA Sequencing
The 150 bp paired-end reads obtained from the RNA-Seq analysis were trimmed using SeqPrep and Sickle to remove adaptor sequences, low-quality reads (those with ambiguous nucleotides and quality scores < 20). The filtered reads were mapped to the genome using HISAT2 v2.0.4. Transcripts were assembled and reconstructed using StringTie based on the HISAT2 mapping files. Several databases were used to annotate gene functions, including the Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), KOG/COG (Eukaryotic Ortholog Groups/Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KO (KEGG Ortholog) and GO (Gene Ontology) databases. Gene expression was normalized using the FPKM (fragments per kilobase of transcript million mapped reads) method. The differentially expressed genes (DEGs) across all the samples were identified using DESeq2 (Love et al., 2014). Genes showing at least a twofold gene expression change with an FDR value < 0.05 were considered significantly differentially expressed. GO (Fisher’s exact test with corrected P-value < 0.05) and KEGG (Fisher’s exact test with corrected P-value < 0.05) enrichment analyses of DEGs associated with the significantly up- and downregulated genes between different stages were performed using Goatools (Klopfenstein et al., 2018) and R packages (Oksanen et al., 2011), respectively. Venn, PCA and hierarchical clustering for all correlation analyses were conducted using R packages. Weighted gene coexpression network analysis (WGCNA) was performed according to the method published previously (Li et al., 2019).
Quantitative Reverse Transcription (RT)-Polymerase Chain Reaction (PCR) and Semiquantitative RT-PCR
For each sample, 1 μg of RNA was used for cDNA synthesis using a TAKARA PrimeScript RT Reagent Kit. Semiquantitative SYBR green-based RT-PCR was performed using SYBR Premix Ex TaqII (Tli RNaseH Plus). Alpha-Tubulin was used as the internal control for normalization. Detailed information for the primers used in the present study is listed in Supplementary Table 1. The presented data are representative of three biological replicates and four technical replicates for each sample. Relative gene expression levels were calculated with the comparative threshold cycle method using a StepOne Plus Real Time PCR System (Applied Biosystems, United States).
Optimization of Culture Medium for Phlebopus portentosus
The medium was optimized based on the results from the genome and transcriptome for suitable carbon sources selection. Stock cultures were grown on complete yeast extract medium (CYM) plates with 2% glucose, 0.2% yeast extracts, 0.2% peptone, 0.1% K2HPO4, 0.05% MgSO4, 0.046% KH2PO4, and 2% agar. All optimization experiments were carried out in 9-cm-diameter Petri dishes (16 ml/petri dish) with modified CYM medium containing 2% different carbon sources. All the samples cultivated at 30°C in the dark for 30 days. The measurement of hyphal diameters was conducted. The significance of growth rate on different carbon sources was determined using t-test.
Results
Genome Assembly
Two sequencing libraries (ONT and Illumina) were constructed for the monokaryotic strain for genome assembly. A total of 951,686 long reads (11,348,781,309 bp) were generated using the ONT platform. The N50 and N90 values of the raw data were 18,704 and 6,226 bp, respectively. After filtering the low-quality reads, 792,451 reads (10,603,105,695 bp) were obtained, with N50, N90 and a mean length of reads 18,866, 7,033 and 13,380 bp. The total data provided 323-fold coverage of the genome (Supplementary Table 2). Illumina sequencing yielded a total of 5.56 Gb clean data, and 20,876,612 clean reads (6,208,784,986 bp) for the monokaryotic strain were obtained by RNA-Seq.
An assembly of 32,742,503 bp across 62 scaffolds with a GC content of 48.92% was obtained (Table 1). Based on the results from the BUSCOO pipeline, only 9 of 290 single-copy entries were missing, suggesting > 95.17% genome completeness. Finally, 9,464 putative genes were predicted using a combination of different pipelines (Table 1). A total of 3393, 3175, 5143, 6751, 5552, 9023, and 9231 encoded proteins were identified with homologous sequences in the GO, KEGG, KOG, Pfam, SwissProt, TrEMBL, and Nr databases, respectively (Table 1). The 103 tRNA genes are dispersed in the genome. However, no rRNA genes were identified. Compared with other fungi used in this study, the genome size of P. portentosus was the smallest and the number of genes was the least except Tuber melanosporum (Table 2).
Table 1. Statistical information associated with the genome sequenced in this study.
Table 2. The total number of CAZyme families in different fungal genomes.
Genome Annotations
Gene annotation and functional categorization were performed using EuKaryotic Orthologous Group (KOG), where 5143 genes were redundantly assigned into 25 categories (Table 1 and Supplementary Figure 1). Except for genes related to general function prediction only (R) and unknown function (S), the number of genes involved in posttranslational modification, protein turnover, and chaperones (O); signal transduction mechanisms (T); and translation, ribosomal structure and biogenesis (J) were the highest, reaching 551, 403, and 334, respectively. The genes assigned to defense mechanisms, nuclear structure, cell motility and extracellular structures were the least abundant, with only 37, 28, 6, and 5 identified, respectively.
In addition, 35.85% of total genes had at least one GO annotation. Within the cellular component (CC) category, 4604 genes were assigned to 11 subcategories, the most abundant being “cell” and “cell part” (Table 1 and Supplementary Figure 2). A total of 3035 genes were classified into 13 molecular function (MF) categories, the greatest number of which was catalytic activity, followed by binding. A total of 1,469 predicted genes were assigned to 13 biological process (BP) GO terms, the most heavily represented being metabolic process and cellular process.
A total of 3,175 genes involved in 108 pathways were detected by KEGG annotation. The most enriched pathways included biosynthesis of amino acids (111), ribosomes (111), RNA transport (109), carbon metabolism (98), spliceosome (92), and purine metabolism (84) (Table 1 and Supplementary Figure 3).
Carbohydrate Metabolism Enzymes Families
A total of 301 CAZyme-coding genes for plant cell wall degradation were detected in the P. portentosus genome. These genes were divided into six families, including 57 auxiliary activities (AA) family genes, 19 carbohydrate-binding module (CBM) family genes, 47 carbohydrate esterase (CE) family genes, 110 glycoside hydrolase (GH) family genes, 62 glycosyl transferase (GT) family genes, and 6 polysaccharide lyase (PL) family genes (Table 2). The total number of CAZyme family genes in P. portentosus was within the range of that observed for mycorrhizal fungi (Table 2). The average number (303) of CAZyme families detected in the genomes of ectomycorrhizal fungi was much lower (ranging from 215 in T. melanosporum to 417 in Piloderma croceum) than that of saprophytic fungi (average number 542, ranging from 449 in A. bisporus to 606 in P. ostreatus) (Table 2). The number of AAs in P. portentosus (57 copies) was similar to that observed in Paxillus involutus (62 copies). P. portentosus had the same number of genes (19) related to CBM families as Amanita muscaria (Table 2). The number of genes belonging to GH families predicted in P. portentosus, A. muscaria, and P. rubicundulus was 110, 109, and 108, respectively. No more than 10 copies of PL genes in mycorrhizal fungi were identified. Except for GT families, the number of AA-, CBM-, CE-, GH-, and PL-encoding genes was much lower in ectomycorrhizal fungi than in saprophytic fungi (Table 2).
Based on the compositions of CAZyme families, the clustering results revealed two clades that included either mycorrhizal or saprophytic fungi. The pattern of CAZyme genes in P. portentosus was similar to that observed in P. rubicundulus, P. involutus, and Suillus luteus (Figure 1). The top five most abundant families in P. portentosus were CE10 (36 copies), GH16 (17 copies), AA7 (16 copies), AA1_1 (12 copies), and CBM5 (8 copies). CE10 genes encode esterases, GH16 enzymes act on xyloglucan and chitin and AA7 genes encode gluco- or chito-oligosaccharide oxidases, all of which were abundant in mycorrhizal and saprophytic fungi (Supplementary Table 1). AA1_1 enzymes are multicopper oxidases and represent the most abundant AA family in P. portentosus and mycorrhizal fungi assayed in the present study, while AA3_2 enzymes belonging to the glucose-methanol-choline (GMC) oxidoreductase family were the most abundant in saprophytic fungi (Supplementary Table 1). In addition, the families GH6 and GH7, which are involved in attacking crystalline cellulose, were present in all saprophytic fungi assayed in the present study and were absent in P. portentosus and other ectomycorrhizal fungi (Supplementary Table 1). The families GH13, GH15, GH27, and GH71 were enriched in P. portentosus, of which the copy numbers of the starch degradation-associated GH13, GH27, and GH71 families were higher than those observed in other mycorrhizal fungi (Supplementary Table 1).
Figure 1. Heatmap clustering of mycorrhizal and saprophytic fungi based on genes related to CAZyme families and lignocellulose decomposition. (A) Heatmap clustering based on CAZyme families; (B) Heatmap clustering based on genes related to lignocellulose decomposition.
Genes Involved in Lignocellulose Decomposition
Low copy numbers of genes were identified belonging to the cellulase, hemicellulase, pectinase, and lignin oxidase families, which directly act on lignocellulose and are listed in Table 3. Hierarchical clustering using the average-linkage method based on the composition of lignocellulases also resulted in the identification of two clades, revealing the different compositions of lignocellulases between ectomycorrhizal and saprophytic fungi (Figure 1B).
Table 3. The distribution of genes related to lignocellulose degradation in P. portentosus and other fungi.
Four cellulase genes were identified in the P. portentosus genome, including 1 endo-beta-1,4-glucanase and 3 β-glucosidases. However, 1,4-β-cellobiosidase, which specifically hydrolyzes 1,4-beta-D-glucosidic linkages, was absent in P. portentosus and other ectomycorrhizal fungi, suggesting that P. portentosus may also have a low capacity to degrade celluloses. A low copy number of genes related to hemicellulases was observed in P. portentosus, with only five hemicellulases (no higher than seven in other ectomycorrhizal fungi) detected compared to saprophytic fungi (ranging from 10 to 29). In particular, no hemicellulases were detected in S. citrinum. α-Glucuronidases were absent in all ectomycorrhizal fungi evaluated in the present study, whereas they were present in all assayed saprophytic fungi. In P. portentosus, no pectinase genes were identified, as was observed in A. muscaria. The number of pectinases in ectomycorrhizal fungi ranged from 0 to 6. In saprophytic fungi, the copy number of pectinases varied significantly in different species, e.g., 3 present in Coprinus cinereus and 16 present in Volvariella volvacea. In addition, pectate lyases, which were detected in all saprophytic fungi, were absent in ectomycorrhizal fungi.
Six multicopper oxidases and one manganese peroxidase were identified in P. portentosus, fewer than was observed in other ectomycorrhizal fungi belonging to Basidiomycota. The number of multicopper oxidases encoded by these fungi was not in agreement with their associated niches and detected in all the fungi. The same pattern of auxiliary enzymes (e.g., aryl-alcohol oxidase, glyoxal oxidase, and benzoquinone reductase) was observed in all the fungi assayed in the present study. Pyranose oxidase, vanillyl-alcohol oxidase and galactose oxidase were not detected in any of the ectomycorrhizal fungi, including P. portentosus, except for L. bicolor.
Gene Expression Profiles in Different Stages
To determine the gene expression patterns during P. portentosus development, 9 samples from 3 growth stages were used for RNA-Seq. The average number of raw reads was 51,592,958, ranging from 47,344,648 in sample S3 to 56,606,312 in sample F3 (Supplementary Table 3). After filtering, the average numbers of clean reads and total bases was 51,161,401 and 7,631,017,892, ranging from 46,876,694 and 6,987,544,020 in sample S3 to 56,203,088 and 8,367,222,043 in sample F3, respectively (Supplementary Table 3). More than 87.62% of the reads within each replicate could be mapped to the reference genome (Supplementary Table 3).
Using an FPKM cutoff value of 1, 8,556, 8,628 and 8,673 genes expressed in stages S, P, and F accounted for 90.40, 91.17, and 91.64% of the total predicted genes, respectively (Supplementary Figure 4A). A total of 8,316 common genes were expressed in all three stages. According to the FPKM values in different stages (FPKM < 10, 10–100, and >100), three gene categories were determined (Supplementary Figure 4B). The majority of genes in all the samples showed moderate expression, with FPKM values ranging from 10 to 100 (Supplementary Figure 4B). The F stage had the greatest number of expressed genes, while the S stage had the fewest. However, the F stage had the fewest highly expressed genes, while the P stage had the most. To determine the correlations among different stages, the nine samples were clustered into three groups based on principal component analysis, which revealed different patterns of gene expression across all stages (Supplementary Figure 5).
Differentially Expressed Genes Across Different Developmental Stages
Based on FPKM values (>1), a corrected P-value of 0.05 and log2 (fold change) of 1 were set as the threshold for significantly different gene expression. A total of 4,921 DEGs were identified in comparisons between the three groups (S vs. P, S vs. F, and P vs. F) (Figure 2). The number of DEGs between the S and P groups was the highest (3600), followed by the S vs. F (2508) and P vs. F (2235) comparisons (Figure 2B). There were 353 common DEGs in all three growth stage comparisons, which were enriched in the carbohydrate metabolism and proteasome categories (Figure 2B). The most unique DEGs (919) were observed between the S and P groups compared to the S vs. F (509) and P vs. F (422) comparisons (Figure 2A). The numbers of downregulated genes in the S vs. P, S vs. F and P vs. F comparisons were higher than those of upregulated genes (Figure 2A).
Figure 2. Analysis of DEGs between different developmental stages. (A) DEG distribution between two analyzed samples. The number of differentially expressed genes is indicated at the top of the histograms. (B) Venn diagrams comparing shared DEGs between the adjacent growth stages.
Functional Classification of Differentially Expressed Genes
All the DEGs were assigned into three different functional GO categories and 41 subcategories, including biological process (BP, 15 subcategories), cellular component (CC, 11 subcategories) and molecular function (MF, 15 subcategories) (Figure 3). The most enriched DEGs were assigned to catalytic activity (GO: 0003824) in all three comparisons. The DEGs in the S vs. P comparison were significantly enriched in drug metabolic process, mitochondrial protein complex and purine-containing compound biosynthetic process. Functional categories for the P vs. F comparison were significantly enriched for gene families associated with cytoplasmic part, organic substance biosynthetic process, cellular biosynthetic process and organic substance biosynthetic process.
Figure 3. GO functional classification of differentially expressed genes. The Y axis represents the GO terms, the upper X axis represents the number of DEGs assigned to GO terms, the basal X axis represents the p-adjusted value, the broken lines represent the number of DEGs in each GO term, and the bars represent the p-adjusted value. (A) The DEGs assigned to GO terms in the S vs. P comparison; (B) indicated the DEGs assigned to GO terms in the P vs. F comparison; (C) the DEGs assigned to GO terms in the S vs. F comparison.
To further understand the molecular and biological functions of the DEGs, they were mapped to the KEGG database (Figure 4). Pathway enrichment analysis revealed that proteasome and DNA replication were the most enriched in the S vs. P comparison. The ribosome pathway was the most significant pathway in the P vs. F comparison, with a P-value of nearly 0, following by DNA replication, oxidative phosphorylation, citrate cycle and peroxisome. Amino sugar and nucleotide sugar metabolism, peroxisome and cysteine and methionine metabolism were the most enriched in the S vs. F comparison.
Figure 4. KEGG pathways of differentially expressed genes. The Y axis represents the KEGG pathways, the upper X axis represents the number of DEGs, the basal X axis represents the p-adjusted value, the broken lines represent the number of DEGs in each KEGG pathway, and the bars represent the p-adjusted value. (A) KEGG pathways of differentially expressed genes in the S vs. P comparison; (B) KEGG pathways of differentially expressed genes in the P vs. F comparison; (C) KEGG pathways of differentially expressed genes in the S vs. F comparison.
Gene Coexpression Network Construction
WGCNA was conducted to uncover the coexpression profiles in successive developmental stages (Figure 5). Eight different modules were identified with high correlation coefficients according to the similarity of expression patterns (Supplementary Tables 5–12). To identify key genes specific to the S stage, the yellow module was analyzed and included a total of 858 genes. Amino sugar and nucleotide sugar metabolism, pyrimidine metabolism, biosynthesis of amino acids, purine metabolism and starch and sucrose metabolism were the most enriched pathways (Supplementary Table 5). Regarding the amino acids biosynthesis pathways, those involved in phenylalanine, tyrosine and tryptophan biosynthesis, histidine metabolism and glycine, serine and threonine metabolism were the most abundant (Supplementary Table 5). The turquoise and red modules (harboring 1,275 and 433 identified genes, respectively) were specific to the P stage, the number of genes in which was higher than was observed in other two stages. In this stage, RNA transport, cell cycle (DNA replication and meiosis), proteasome and spliceosome were the most enriched pathways (Supplementary Tables 6,7). The blue (63 genes) and black (63 genes) modules were highly associated with the F stage. Biosynthesis of amino acids, glutathione metabolism, peroxisome, glycerophospholipid metabolism and MAPK signaling pathway were the most enriched in this stage (Supplementary Tables 8,9).
Figure 5. WGCNA of genes in different developmental stages in P. portentosus. Module-trait weighted correlations and corresponding P-values for the identified gene module and their developmental stages. The label color on the right represents the strength of correlation, from 1 (red) to –1 (blue). Each column corresponds to a developmental stage. The eight modules and the number of module member genes are shown on the left panel. The boxes with different color indicated the significant correlation between the module eigengene and the stage.
Differential Expression of Carbohydrate Metabolism Enzymes Families and Lignocellulose Decomposition-Related Genes
Of the 301 CAZyme genes predicted from the genome, 277, 286 and 288 genes were expressed in the S, P, and F stages, respectively (Supplementary Table 4). The most highly expressed CAZyme genes were different in the three stages (Supplementary Table 4). In the S stage, the top five families with the highest FPKM values were CBM13 (1872.75), GH16 (1230.55), GH128 (1168.78), CBM13 (923.36) and AA5_1 (532.96). In the P stage, AA2 (1391.67), GT4 (992.50), GH128 (867.38), AA3_2 (298.65) and CBM50 (288.51) were more highly expressed than other families. The families GH16 (6346.44), GH18 (1097.48), GH128 (1003.62), GT4 (969.11), and GT2_Glyco_trans_2_3 (660.98) were overrepresented in stage F. The high expression of the glycoside hydrolase families GH16, GH18 and GH128 in all three stages indicated that P. portentosus has a high potential for starch or xylan degradation.
A total of 106 and 78 DEGs annotated as CAZyme families were detected in the S vs. P and P vs. F comparisons, respectively (Supplementary Table 4). In the S vs. P comparison, out of the 106 identified DEGs, 46 were upregulated and 60 were downregulated. Compared to the S stage, CBM13, GH16, AA5_1, GH18, AA6, AA3_2, GH17, GH5_30 and GT2_Chitin_synth_1 were significantly downregulated and AA2, GT4 and GH71 were upregulated in the P stage. In the P vs. F comparison, 52 DEGs were upregulated, including GH16, GH18, GT2_Glyco_trans_2_3, GH13_1 and GH31, while 26 DEGs were downregulated, including AA2, GH16, CE10, GH16, and CE4.
Differential Expression of Genes Involved in Lignocellulose Decomposition
Twenty-four genes involved in lignocellulose decomposition were predicted from the P. portentosus reference genome, of which 15 genes were differentially expressed in the different developmental stages (Table 4). The most highly expressed genes involved in lignocellulose decomposition in all 3 stages included lignin oxidases (manganese peroxidase) and lignin-degrading auxiliary enzymes (alcohol oxidase, benzoquinone reductase), suggesting that P. portentosus has the potential to degrade lignin (Table 4). However, endo-beta-1,4-glucanase, which acts on cellulose, was expressed at low levels. Compared to that observed in the S stage, the expression of manganese peroxidase (EVM0004344) in the P stage was significantly upregulated (293.84 vs. 1391.67, respectively, a nearly fivefold increase) (Table 4), the expression of which subsequently decreased in the F stage. The RT-PCR results confirmed the observed expression profiles (Figure 6A). Multicopper oxidase laccases (EVM0002770, EVM0008799, and EVM0004198) were expressed at the highest levels in the S stage (Table 4 and Figure 6B). The lignin-degrading auxiliary enzymes glyoxal oxidase (EVM0006667), benzoquinone reductase (EVM0006980 and EVM0009178) and glucose oxidase (EVM0003870) were more highly expressed in the P stage than in the S and F stages (Table 4). The expression of aryl-alcohol oxidase (EVM0008750) gradually increased with the development of fruiting bodies, peaking in the F stage (Table 4). Endo-beta-1,4-glucanase expression was lower in the S (11.68) and F (23.71) stages than in the P stage (46.76) (Table 4). The expression of β-glucosidase (EVM0002244) increased from stages P (29.79) to F (105.70) (Table 4 and Figure 6D). Endo-1,4-betaxylanase (EVM0003710), which acts on hemicellulose, was expressed at higher levels in the F stage than in the P and S stages (Table 4 and Figure 6C).
Table 4. FPKM of DEGs annotated as lignocellulose degradation enzymes.
Figure 6. Relative expression of some lignin oxidases and cellulases in different developmental stages of P. portentosus using RT-PCR. (A) Manganese peroxidase; (B) multicopper oxidase; (C) endo-beta-1,4-glucanase; (D) β-glucosidase. Different lowercase letters above bars indicate a significant difference according to the Duncan’s multiple range test (p < 0.05).
Optimization of Culture Medium
From the genome and transcriptome results, P. portentosus exhibited a limited ability to degrade lignocellulose. The gene families acting on starch, xylan and chitin were most abundant in the genome and showed high expression in all the stages. Unique carbon sources (starch, xylan, and chitin) were further used to screen the most suitable culture medium, glucose as a control. By measuring the growth diameter of mycelia after 30 days of cultivation, the growth rate of P. portentosus on starch and xylan was much faster than that on glucose and chitin (Figure 7A). In addition to the growth rate, the morphology of the colonies also differed. Aerial mycelia of P. portentosus were much thicker, and pigmentation was stronger on glucose and starch than on xylan (Figure 7B). However, the integrity of mycelial edges were much better on xylan than on glucose and starch (Figure 7B). Thus, xylan and starch may be more suitable carbon sources for P. portentosus than chitin.
Figure 7. The growth rate and morphology of P. portentosus cultured on different carbon sources. (A) Growth rate; (B) morphology of P. portentosus on different carbon sources after 30 days of cultivation, (B1) glucose, (B2) starch, (B3) xylan, (B4) chitin. Different lowercase letters above bars indicate a significant difference according to the Duncan’s multiple range test (p < 0.05).
Discussion
As the first ectomycorrhizal fungus in the order Boletaceae to have been industrially cultivated at a large scale, P. portentosus has attracted extensive attention worldwide (Zhang et al., 2017). However, the molecular mechanisms associated with the ability of P. portentosus to promote plant cell wall degradation and fruiting body development remain ambiguous. Genomics and transcriptomics provide an unprecedented means of elucidating the developmental and metabolic biological aspects of mushrooms. In the present study, genome and transcriptome sequencing was conducted to identify CAZyme family genes and DEGs in different developmental stages of P. portentosus. The results revealed the presence of no or few copies of carbohydrate metabolism enzymes (CAZymes) acting on cellulose, xylan, pectin and lignin in P. portentosus, an ectomycorrhizal fungus (Kohler et al., 2015). Some conserved pathways related to the proteasome, DNA replication, ribosome pathway and amino acid metabolism processes were correlated with the development of fruiting bodies.
Ectomycorrhizal fungi are considered mutualists that trade host photoassimilates for nutrients (Clavel et al., 2021). ECM fungi have been reported to have evolved from saprotrophic fungi and exhibit an extensive loss of enzymatic repertoires for lignocellulose decomposition during the transition to mycorrhizal habits (Kohler et al., 2015). The use of coculture and isotopic methods revealed that P. portentosus can form ectomycorrhizal symbioses with plants (Kumla et al., 2016). However, this fungus can be artificially cultivated in vitro, and an improved understanding of the repertoires of PCWDEs in P. portentosus is needed to elucidate the mechanisms associated with nutrient acquisition in vitro.
The composition of lignocellulose-degrading systems harbored by P. portentosus revealed ECM characteristics. Members of the families GH6 and GH7, which degrade crystalline cellulose (Baldrian and Valáková, 2008; Floudas et al., 2012), are absent in P. portentosus and other ECM fungi assayed in the present study except Piloderma croceum (1 GH7 copy). In addition to GH6 and GH7, fewer AA9 lytic polysaccharide monooxygenases (2–5 copies, LPMOs) were present in all the genomes of ectomycorrhizal fungi assayed in the present study than in white-rot or straw-rot fungi (11–34). These results are in agreement with reports of ECM fungi belonging to Boletales, which have lost nearly all copies of GH6 and GH7 genes and have no more than five copies of AA9 LPMOs (Kohler et al., 2015). These results also suggested that ECM fungi, including P. portentosus, have a limited capacity to degrade cellulose. In addition, no genes related to pectin degradation were detected in P. portentosus and A. muscaria. The lack of pectinases in ECM fungi has also been demonstrated using omics analyses (Kohler et al., 2015). Genes encoding enzymes mediating the degradation of pectin were shown to be reduced or completely lost among EM Amanita (Jaqueline et al., 2018). Each species used in the study had a different repertoire of multicopper oxidase (MCO) genes, ranging from 4 to 13 in saprophytic fungi and from 5 to 19 in ECM fungi (Bödeker et al., 2009; Kohler et al., 2015; Shah et al., 2016). Six laccases (MCO) were identified in the P. portentosus genome, higher than that observed for S. commune but lower than other ECM fungi (except T. melanosporum). In addition to the low copy number of lignocellulose-degrading genes observed in the P. portentosus genome, the RNA-seq analysis of the different stages showed the low expression of these genes, with measured FPKM values of half of these lignocellulose-degrading genes less than 10. These results were also in agreement with the characteristics of ectomycorrhizal fungi revealed by omics analyses in other studies, and ectomycorrhizal fungi have low potential to degrade lignocellulose (Martin et al., 2016; Shah et al., 2016; Miyauchi et al., 2020).
Although P. portentosus was shown to have few copies and low expression of lignocellulose enzymes, the expression of these genes was also dependent on the stage of mushroom development. The expression of endo-beta-1,4-glucanase, β-glucosidase and endo-1,4-betaxylanase was highest in the P or F stages, which may suggest that cellulase and hemicellulase genes are potentially correlated with fruiting body formation and maturation. These results are in agreement with those obtained for some white-rot or straw-rot mushrooms, e.g., A. bisporus, P. ostreatus, L. edodes, L. tigrinus, and G. lucidum (Shoji and Royse, 2001; Lechner and Papinutti, 2006; Elisashvili et al., 2008; Song et al., 2018; Wu et al., 2018; Zhou et al., 2018). The activities of multicopper oxidases (laccases) and most lignin-degrading auxiliary enzymes reached their maximum levels in the M stage, similar to other mushrooms (Shoji and Royse, 2001; Lechner and Papinutti, 2006; Elisashvili et al., 2008), suggesting the potential to degrade lignin during the mycelial growth phase. Manganese peroxidase activity was reported to be high during the colonization stage and decrease during the first primordia and fruiting body formation stages in some saprophytic mushrooms (Velázquez-Cedeño et al., 2002; Elisashvili et al., 2008). However, the relative expression of manganese peroxidases in P. portentosus was inconsistent with previous reports, reaching the highest levels in the primordia stage and then decreasing in the fruiting body formation stage. The expression patterns of these enzymes in P. portentosus need to be further studied.
Some conserved and specific pathways have been shown to be enriched in the different stages of mushroom development (Wang et al., 2013, 2018; Zhang et al., 2015; Song et al., 2018; Xie et al., 2018; Zhou et al., 2018; Hao et al., 2019; Wei et al., 2019; Tang et al., 2020). Starch and sucrose metabolism pathways were enriched in the mycelial stage in P. portentosus, which were also specific to the hyphal stage of Chinese Cordyceps (Li et al., 2019). Starch and sucrose may be important carbon sources for the growth of mycelia. Cultivation using starch and xylan as substrates demonstrated that metabolism plays important roles in acquiring nutrient acquisition in P. portentosus. Some ectomycorrhizal fungi belonged to Amanita, Boletus, Laccaria, and Tricholoma revealed the ability to utilize starch as a growth substrate (Ohta, 1997; Kusuda et al., 2006). Tuber maculatum mycelium utilized xylan as energy source during the growth cycle (Bedade et al., 2017). Both DNA replication and proteasome were enriched in P and F stages in P. portentosus, which were detected to be associated with development and maturation in other mushrooms (Masato et al., 2010; Rodenburg et al., 2018; Song et al., 2018; Wang et al., 2018). DNA replication is important for premeiotic replication, karyogamy, and meiosis during maturation in mushrooms, e.g., L. edodes (Song et al., 2018; Wang et al., 2018) and Botrytis cinerea (Rodenburg et al., 2018). The proteasome is presumably linked with fruiting body development in basidiomycetes by degrading unnecessary proteins through the ubiquitin-proteasome pathway (Masato et al., 2010). Furthermore, peroxisome, ribosome, oxidative phosphorylation, citrate cycle and MAPK signaling pathways were enriched in the stage of F in P. portentosus. Ribosome pathways were increased significantly during the fruiting body growth stage in Ophiocordyceps sinensis (Zhang et al., 2021), F. velutipes (Liu et al., 2014), Auricularia polytricha (Zhou et al., 2014), and Agaricus blazei (Lu et al., 2020) to provide proteins for fruiting body formation or carry out some other functions, e.g., DNA repair, development, and cell division. In fungi, the activity of peroxisomes was required for sexual development (Peraza-Reyes and Berteaux-Lecellier, 2013). Oxidative phosphorylation and citrate cycle were the key metabolic pathways for the growth and development of Clavariadelphus pistillaris (Tang et al., 2019). Higher expression of MAPK genes involved in MAPK signaling pathways is required for fruiting in some fungi, which plays an important role in promoting sexual reproduction, regulating the polar growth of mycelia and rebuilding of the cell wall in the mature fruiting body (Pöggeler et al., 2006; Zhang et al., 2015; Li et al., 2019).
Conclusion
Our results revealed fewer cell wall-degrading enzymes encoded in the genome of P. portentosus than that observed in some saprophytic mushrooms with specific ectomycorrhizal niches. The low expression of plant cell wall-degrading genes confirmed the limited ability of P. portentosus to degrade lignocellulose. In addition, some conserved and specific pathways were enriched in the different developmental stages of P. portentosus, e.g., starch and sucrose metabolism enriched in the mycelial stage and DNA replication, proteasome and MAPK signaling pathways enriched in the fruiting body stage. The results of the present study advance our understanding of the molecular mechanisms of P. portentosus lignocellulose degradation and fruiting body development.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA743691 (SRR15039576–SRR15039584). https://www.ncbi.nlm.nih.gov/, JAHRGP000000000.
Author Contributions
R-HY, D-PB, K-PJ, S-ZL, and QT contributed to the conception and design of the study. J-NW, YL, S-ZL, YC, and G-YJ prepared the samples. R-HY and K-PJ performed the analysis and wrote the manuscript draft. All authors approved the final version of the manuscript.
Funding
This work was supported by Shanghai Action Plan for Scientific and Technological Innovation (19391901400) and the National Natural Science Foundation of China (32072644).
Conflict of Interest
K-PJ, YC, S-ZL, and G-YJ are employed by Hongzhen Agricultural Science and Technology Co. Ltd.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2021.732458/full#supplementary-material
Footnotes
References
Baldrian, P., and Valáková, V. (2008). Degradation of cellulose by basidiomycetous fungi. FEMS Microbiol. Rev. 32, 501–521.
Bedade, D., Berezina, O., Singhal, R., Deska, J., and Shamekh, S. (2017). Extracellular xylanase production from a new xylanase producer Tuber maculatum mycelium under submerged fermentation and its characterization. Biocatalysis Agric. Biotechnol. 11, 288–293.
Blanco, E., Parra, G., and Guigo, R. (2002). “Using GeneID to identify genes,” in Current Protocols in Bioinformatics, ed. A. Baxevanis (New York, NY: John Wiley & Sons Inc).
Bödeker, I. T. M., Nygren, C. M. R., Taylor, A. F. S., Olson, A., and Lindahl, B. D. (2009). ClassII peroxidase-encoding genes are present in a phylogenetically wide range of ectomycorrhizal fungi. ISME J. 3, 1387–1395. doi: 10.1038/ismej.2009.77
Boeckmann, B., Blatter, M. C., Famiglietti, L., Hinz, U., Lane, L., Roechert, B., et al. (2005). Protein variety and functional diversity: swiss-prot annotation in its biological context. Comptes Rendus Biol. 328, 882–899.
Burge, C. B., and Karlin, S. (1998) Finding the genes in genomic DNA. Curr. Opin. Struct. Biol. 8, 346–354.
Cao, Y., Zhang, Y., Yu, Z., Wang, P., Tang, X., He, X., et al. (2015). Genome sequence of Phlebopus portentosus strain PP33, a cultivated Bolete. Genome Announcements 3, e326–e315.
Chang, S. T., Li, G., and Peberdy, J. F. (1985). Isolation of protoplasts from edible fungi. Mircen J. Appl. Microbiol. Biotechnol. 1, 185–193.
Chen, L., Gong, Y., Cai, Y., Wei, L., and Bian, Y. (2016). Genome sequence of the edible cultivated mushroom Lentinula edodes (shiitake) reveals insights into lignocellulose degradation. PLoS One 11:e0160336. doi: 10.1371/journal.pone.0160336
Clavel, J., Lembrechts, J., Alexander, J., Haider, S., Lenoir, J., Milbau, A., et al. (2021). The role of arbuscular mycorrhizal fungi in nonnative plant invasion along mountain roads. New Phytol. 230, 1156–1168. doi: 10.1111/nph.16954
Conesa, A., Gotz, S., Garcia-Gomez, J. M., Terol, J., Talon, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610
Elisashvili, V., Kachlishvili, E., and Penninckx, M. (2008). Lignocellulolytic enzymes profile during growth and fruiting of Pleurotus ostreatus on wheat straw and tree leaves. Acta Microbiol. Immunol. Hungarica 55, 157–168. doi: 10.1556/AMicr.55.2008.2.7
Finn, R. D., Mistry, J., Tate, J., Coggill, P., Heger, A., Pollington, J. E., et al. (2010). Pfam protein families database. Nuclc Acids Res. 40, 290–301.
Floudas, D., Binder, M., Riley, R., Barry, K., Blanchette, R. A., Henrissat, B., et al. (2012). The paleozoic origin of enzymatic lignin decomposition reconstructed from 31 fungal genomes. Science 336, 1715–1719. doi: 10.1126/science.1221748
Haas, B. J., Salzberg, S. L., Zhu, W., Pertea, M., Allen, J. E., Orvis, J., et al. (2008). Automated eukaryotic gene structure annotation using EVidenceModeler and the Program to Assemble Spliced Alignments. Genome Biol. 9:R7. doi: 10.1186/gb-2008-9-1-r7
Hall, I. R., Yun, W., and Amicucci, A. (2003). Cultivation of edible ectomycorrhizal mushrooms. Trends Biotech. 21, 433–438.
Hammel, K. E. (1997). “Fungal degradation of lignin,” in Driven by Nature: Plant Litter Quality and Decomposition, eds G. Cadisch and K. E. Giller (Wallingford: CAB International), 33–45.
Hao, H., Zhang, J., Wang, H., Wang, Q., and Chen, H. (2019). Comparative transcriptome analysis reveals potential fruiting body formation mechanisms in Morchella importuna. AMB Express 9:103. doi: 10.1186/s13568-019-0831-4
Huang, L., Zhang, H., Wu, P., Sarah, E., Li, X., Tanner, Y., et al. (2018). dbCAN-seq: a database of carbohydrate-active enzyme (CAZyme) sequence and annotation. Nucleic Acids Res. 46, D516–D521.
Isikhuemhen, O. S., and Mikiashvilli, N. A. (2009). Lignocellulolytic enzyme activity, substrate utilization, and mushroom yield by Pleurotus ostreatus cultivated on substrate containing anaerobic digester solids. J. Industrial Microbiol. 36, 1353–1362. doi: 10.1007/s10295-009-0620-1
Jaqueline, H., Inger, S., Maryam, C. D. M., Matthieu, H., Bernard, H., and Anne, P. (2018). Rapid divergence of genome architectures following the origin of an ectomycorrhizal symbiosis in the genus Amanita. Mol. Biol. Evol. 35, 2786–2804. doi: 10.1093/molbev/msy179
Jens, K., Frank, H., and Jan, G. (2019). GeMoMa: homology-based gene prediction utilizing intron position conservation and RNA-Seq data. Methods Mol. Biol. 1962, 161–177.
Ji, K. P., Cao, Y., Zhang, C. X., He, M. X., Liu, J., Wang, W. B., et al. (2011). Cultivation of Phlebopus portentosus in southern China. Mycol. Progr. 10, 293–300.
Kim, D., Paggi, J. M., Park, C., Bennett, C., and Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37:1. doi: 10.1038/s41587-019-0201-4
Klopfenstein, D. V., Zhang, L., Pedersen, B. S., Fidel, R., Alex, W. V., Aurélien, N., et al. (2018). GOATOOLS: a python library for gene ontology analyses. Sci. Rep. 8:10872. doi: 10.1038/s41598-018-28948-z
Kohler, A., Kuo, A., Nagy, L. G., Morin, E., Barry, K. W., Buscot, F., et al. (2015). Convergent losses of decay mechanisms and rapid turnover of symbiosis genes in mycorrhizal mutualists. Nat. Genet. 47, 410–415. doi: 10.1038/ng.3223
Kumla, J., Danell, E., and Lumyong, S. (2015). Improvement of yield for a tropical black bolete, Phlebopus portentosus, cultivation in northern Thailand. Mycoscience 56, 114–117.
Kumla, J., Hobbie, E. A., Suwannarach, N., and Lumyong, S. (2016). The ectomycorrhizal status of a tropical black bolete, Phlebopus portentosus, assessed using mycorrhizal synthesis and isotopic analysis. . Mycorrhiza 26, 333–343. doi: 10.1007/s00572-015-0672-1
Kumla, J., Suwannarach, N., Tanruean, K., and Lumyong, S. (2021). Comparative evaluation of chemical composition, phenolic compounds, and antioxidant and antimicrobial activities of tropical black bolete mushroom using different preservation methods. Foods 10:781. doi: 10.3390/foods10040781
Kusuda, M., Ueda, M., Konishi, Y., Araki, Y., Yamanaka, K., Nakazawa, M., et al. (2006). Detection of β-glucosidase as a saprotrophic ability from an ectomycorrhizal mushroom, Tricholoma matsutake. Mycoscience 47, 184–189.
Lechner, B. E., and Papinutti, V. L. (2006). Production of lignocellulosic enzymes during growth and fruiting of the edible fungus Lentinus tigrinus on wheat straw. Process Biochem. 41, 594–598.
Li, X., Wang, F., Liu, Q., Li, Q., Qian, Z., Zhang, X., et al. (2019). Developmental transcriptomics of Chinese cordyceps reveals gene regulatory network and expression profiles of sexual development-related genes. BMC Genomics 20:337. doi: 10.1186/s12864-019-5708-z
Liu, F., Wang, W., and Xie, B. (2014). Comparison of gene expression patterns in the mycelium and primordia of Flammulina velutipes, strain 1123. Acta Edulis Fungi 21, 1–7.
Liu, X. B., Orcid, I., Xia, E. H., Li, M., Cui, Y. Y., Wang, P. M., et al. (2020). Transcriptome data reveal conserved patterns of fruiting body development and response to heat stress in the mushroom-forming fungus Flammulina filiformis. PLoS One 15:e0239890. doi: 10.1371/journal.pone.0239890
Love, M. I, Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8
Lowe, T. M., and Eddy, S. R. (1997). tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 25, 955–964. doi: 10.1093/nar/25.5.955
Lu, Y. P., Liao, J. H., Guo, Z. J., Cai, Z. X., and Chen, M. Y. (2020). Genome survey and transcriptome analysis on mycelia and primordia of Agaricus blazei. BioMed. Res. Int. 11, 1–12. doi: 10.1155/2020/1824183
Majoros, W., Pertea, M., and Salzberg, S. (2004). TigrScan and GlimmerHMM: two open source ab initio eukaryotic gene-finders. Bioinformatics 20, 2878–2879. doi: 10.1093/bioinformatics/bth315
Mario, S., and Burkhard, M. (2005). AUGUSTUS: a web server for gene prediction in eukaryotes that allows user-defined constraints. Nucleic Acids Res. 33, 465–467. doi: 10.1093/nar/gki458
Martin, F., Kohler, A., Murat, C., Veneault-Fourrey, C., and Hibbett, D. S. (2016). Unearthing the roots of ectomycorrhizal symbioses. Nat. Rev. Microbiol. 14, 760–773. doi: 10.1038/nrmicro.2016.149
Masato, Y., Sou, S., Kou, S., Goro, T., Satoshi, I., Mitsuo, O., et al. (2010). Isolation and analysis of genes specifically expressed during fruiting body development in the basidiomycete Flammulina velutipes by fluorescence differential display. FEMS Microbiol. Lett. 254, 165–172. doi: 10.1111/j.1574-6968.2005.00023.x
Miyauchi, S., Kiss, E., Kuo, A., Drula, E., Kohler, A., Sánchez-García, M., et al. (2020). Large-scale genome sequencing of mycorrhizal fungi provides insights into the early evolution of symbiotic traits. Nat. Commun. 11:5125. doi: 10.1038/s41467-020-18795-w
Nawrocki, E. P. (2014). Annotating functional RNAs in genomes using Infernal. Methods Mol. Biol. 1097, 163–197.
Ohta, A. (1997). Ability of ectomycorrhizal fungi to utilize starch and related substrates. Mycoscience 38, 403–408.
Peraza-Reyes, L., and Berteaux-Lecellier, V. (2013). Peroxisomes and sexual development in fungi. Front. Physiol. 4:244. doi: 10.3389/fphys.2013.00244
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Pham, N., Yamada, A., Shimizu, K., Noda, K., Dang, L., and Suzuki, A. (2012). A sheathing mycorrhiza between the tropical bolete Phlebopus spongiosus and Citrus maxima. Mycoscience 53, 347–353.
Pöggeler, S., Nowrousian, M., Ringelberg, C., Loros, J. J., Dunlap, J. C., and Kück, U. (2006). Microarray and real-time PCR analyses reveal mating type-dependent gene expression in a homothallic fungus. Mol. Genet. Genomics 275, 492–503. doi: 10.1007/s00438-006-0107-y
Quevillon, E., Silventoinen, V., Pillai, S., Harte, N., Mulder, N., Apweiler, R., et al. (2005). InterProScan: protein domains identifier. Nucleic Acids Res. 33, W116–W120.
Rédei, G. (2008). TREMBL. Encyclopedia of Genetics, Genomics, Proteomics and Informatics. Berlin: Springer.
Rodenburg, S., Terhem, R. B., Veloso, J., Stassen, J., Kan, J. V., and Turgeon, B. G. (2018). Functional analysis of mating type genes and transcriptome analysis during fruiting body development of Botrytis cinerea. Mbio 9, e1939–e1917. doi: 10.1128/mBio.01939-17
Ruan, J., and Li, H. (2020). Fast and accurate long-read assembly with wtdbg2. Nat. Methods 17, 1–4. doi: 10.1038/s41592-019-0669-3
Sahoo, N. (2017). Sequence Base-Calling Through Albacore Software: A Part of the Oxford Nanopore Technology. MSc thesis. Odisha: Orissa University of Agriculture and Technology, Bhubaneswar, India. Vol. 27, 722–736.
Sanmee, R., Lumyong, P., Dell, B., and Lumyong, S. (2010). In vitro cultivation and fruit body formation of the black bolete, Phlebopus portentosus, a popular edible ectomycorrhizal fungus in Thailand. Mycoence 51, 15–22.
Sergey, K., Brian, P., and Konstantin, B. (2017). Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Research 27, 722–736. doi: 10.1101/gr.215087.116
Shah, F., Nicolás, C., Bentzer, J., Ellström, M., Smits, M., Rineau, F., et al. (2016). Ectomycorrhizal fungi decompose soil organic matter using oxidative mechanisms adapted from saprotrophic ancestors. N. Phytol. 209, 1705–1719. doi: 10.1111/nph.13722
Shoji, O., and Royse, D. J. (2001). Transcriptional regulation of laccase and cellulase genes during growth and fruiting of Lentinula edodes on supplemented sawdust. FEMS Microbiol. Lett. 201, 111–115. doi: 10.1111/j.1574-6968.2001.tb10741.x
Simão, F., Waterhouse, R. M., Panagiotis, I., Kriventseva, E. V., and Zdobnov, E. M. (2015). BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210– 3212.
Song, H. Y., Dae-Hyuk, K., and Jung-Mi, K. (2018). Comparative transcriptome analysis of dikaryotic mycelia and mature fruiting bodies in the edible mushroom Lentinula edodes. Sci. Rep. 8:8983. doi: 10.1038/s41598-018-27318-z
Sun, Z., Hu, M., Sun, Z., Zhu, N., Yang, J., Ma, G., et al. (2018). Pyrrole alkaloids from the edible mushroom Phlebopus portentosus with their bioactive activities. Molecules 23:1198. doi: 10.3390/molecules23051198
Tang, X., Ding, X., Dong, M., Zhu, M., Song, Z., and Hou, Y. (2019). Transcriptome analysis of Clavariadelphus pistillaris fruiting bodies atdifferent development stages. Biotechnol. Bull. 35, 119–129.
Tang, X., Ding, X., and Hou, Y. L. (2020). Comparative analysis of transcriptomes revealed the molecular mechanism of development of Tricholoma matsutake at different stages of fruiting bodies. Food Sci. Biotechnol. 29, 939–951. doi: 10.1007/s10068-020-00732-8
Tseng, M. C., and Luong, J. H. (1984). Mushroom cultivation - technology for commercial production. Annu. Rep. Fermentation Processes 7, 45–79.
Velázquez-Cedeño, M., Mata, G., and Savoie, J. M. (2002). Waste-reducing cultivation of Pleurotus ostreatus and Pleurotus pulmonarius on coffee pulp: changes in the production of some lignocellulolytic enzymes. World J. Microbiol. Biotechnol. 18, 201–207.
Walker, B. J., Abeel, T., Shea, T., Priest, M., and Earl, A. M. (2014). Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS One 9:e112963. doi: 10.1371/journal.pone.0112963
Wang, M., Gu, B., Huang, J., Jiang, S., Chen, Y., Yin, Y., et al. (2013). Transcriptome and proteome exploration to provide a resource for the study of Agrocybe aegerita. PLoS One 8:e56686. doi: 10.1371/journal.pone.0056686
Wang, Y., Zeng, X., and Liu, W. (2018). De novo transcriptomic analysis during Lentinula edodes fruiting body growth. Gene 641, 326–334. doi: 10.1016/j.gene.2017.10.061
Wei, L., Cai, Y., He, P., Chen, L., and Bian, Y. (2019). Comparative transcriptomics reveals potential genes involved in the vegetative growth of Morchella importuna. 3 Biotech 9:81. doi: 10.1007/s13205-019-1614-y
Wu, B., Xu, Z., Knudson, A., Carlson, A., Chen, N., Kovaka, S., et al. (2018). Genomics and development of Lentinus tigrinus: a white-rot wood-decaying mushroom with dimorphic fruiting bodies. Genome Biol. Evol. 10, 3250–3261. doi: 10.1093/gbe/evy246
Xie, C., Gong, W., Zhu, Z., Yan, L., Hu, Z., and Peng, Y. (2018). Comparative transcriptomics of Pleurotus eryngii reveals blue-light regulation of carbohydrate-active enzymes (CAZymes) expression at primordium differentiated into fruiting body stage. Genomics 110, 201–209. doi: 10.1016/j.ygeno.2017.09.012
Yang, R. H., Bao, D. P., Guo, T., Li, Y., and Tan, Q. (2019). Bacterial profiling and dynamic succession analysis of Phlebopus portentosus casing soil using MiSeq sequencing. Front. Microbiol. 10:1927. doi: 10.3389/fmicb.2019.01927
Yang, R. H., Li, Y., Wáng, Y., Wan, J. N., Zhou, C. L., Wãng, Y., et al. (2016). The genome of Pleurotus eryngii provides insights into the mechanisms of wood decay. J. Biotechnol. 239, 65–67. doi: 10.1016/j.jbiotec.2016.10.007
Yuki, M., Masumi, I., Shujiro, O., Yoshizawa, A. C., and Minoru, K. (2007). KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 35, W182–W185.
Zhang, C., He, M., Jing, L., Xu, X., and Yun, W. (2017). Brief introduction to a unique edible Bolete-Phlebopus portentosus in southern China. J. Agric. Sci. Technol. B 7, 386–394.
Zhang, H., Yue, P., Tong, X., Bai, J., and Guo, J. (2021). mRNA-seq and miRNA-seq profiling analyses reveal molecular mechanisms regulating induction of fruiting body in ophiocordyceps sinensis. Sci. Rep. 11:12944. doi: 10.1038/s41598-021-91718-x
Zhang, J., Ren, A., Chen, H., Zhao, M., Shi, L., Chen, M., et al. (2015). Transcriptome analysis and its application in identifying genes associated with fruiting body development in Basidiomycete Hypsizygus marmoreus. PLoS One 10:e0123025. doi: 10.1371/journal.pone.0123025
Zhou, S., Zhang, J., Ma, F., Tang, C., Tang, Q., and Zhang, X. (2018). Investigation of lignocellulolytic enzymes during different growth phases of Ganoderma lucidum strain G0119 using genomic, transcriptomic and secretomic analyses. PLoS One 13:e0198404. doi: 10.1371/journal.pone.0198404
Keywords: ectomycorrhizal fungi, genome, transcriptome, CAZymes, plant cell wall degradation
Citation: Wan J-N, Li Y, Guo T, Ji G-Y, Luo S-Z, Ji K-P, Cao Y, Tan Q, Bao D-P and Yang R-H (2021) Whole-Genome and Transcriptome Sequencing of Phlebopus portentosus Reveals Its Associated Ectomycorrhizal Niche and Conserved Pathways Involved in Fruiting Body Development. Front. Microbiol. 12:732458. doi: 10.3389/fmicb.2021.732458
Received: 29 June 2021; Accepted: 03 September 2021;
Published: 29 September 2021.
Edited by:
Baolei Jia, Chung-Ang University, South KoreaReviewed by:
Shireesh Srivastava, International Centre for Genetic Engineering and Biotechnology, IndiaLujiang Hao, Qilu University of Technology, China
Yee Shin Tan, University of Malaya, Malaysia
Copyright © 2021 Wan, Li, Guo, Ji, Luo, Ji, Cao, Tan, Bao and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kai-Ping Ji, jkpcnlc@126.com; Rui-Heng Yang, yangrh@saas.sh.cn
†These authors have contributed equally to this work