The P. tuber-regium strain ACCC 50657 was kept in the Laboratory of Food Microbiology, Huazhong Agricultural University. Monokaryotic strains 2 and 18 were a pair of compatible single-spore isolates from the fruiting body of P. tuber-regium strain ACCC 50657, and dikaryotic strain 218 was obtained by pairing the two monokaryotic strains. These strains were maintained and subcultured on solid potato dextrose agar (PDA) medium (potato 200 g/L, glucose 20 g/L, and agar 20 g/L) at room temperature and 32°C, respectively.

To induce formation of sclerotia, culture bags were prepared that were composed of cottonseed shell, bran, gypsum, lime, and H₂O. These 300-g sterilized culture bags were inoculated with each P. tuber-regium strain and incubated at 32°C for 15 d to ensure complete colonization. Then, each culture bag was buried in 2.5 L of nutrient soil using a plastic vase (approximately 5.5 L) with a rounded bottom and incubated at room temperature with proper humidity.

Monokaryotic strain 18 was selected for genome sequencing. It was grown on PDA medium, covered with cellophane, and incubated at 32°C for 7 d in darkness. Genomic DNA was extracted using the CTAB extraction protocol (Cota-Sánchez et al., 2006). The DNA quality and quantity were determined using a Qubit dsDNA HS Assay Kit (Invitrogen, United States) and an Agilent Bioanalyzer 2100 (Agilent Technologies, United States).

After filtering low-quality reads, the genome survey was performed by counting the frequency of 17-mers from 5.7 Gb of data. Ten thousand randomly picked read pairs were blasted to the NCBI non-redundant nucleotide (nt) database to check for obvious sample contamination. Subsequently, 10 μg genomic DNA was used for 20 kb template library preparation using the BluePippin Size Selection system (Sage Science, United States) according to the manufacturer’s protocol (Pacific Biosciences, United States). The library was sequenced on the Pacific Biosciences Sequel platform. Assembly was performed using the Canu software with default parameters (Koren et al., 2017), and then, Pilon software was used to correct the preassembled structure according to NGS reads to obtain the final assembly result (Walker et al., 2014).

The protein-coding genes were annotated using predicted proteins to compare with seven major databases (NR, Swiss-Prot, KOG, KEGG, GO, Pfam, and TrEMBL), and the functional prediction of proteins was performed to combine the functional information of the proteins in the databases. Genes with putative CAZymes were annotated by BLASTP analysis according to the Carbohydrate Active Enzymes database¹ (Cantarel et al., 2009).

There are two types of repetitive elements in the genome: tandem (e.g., microsatellites and minisatellites) and interspersed repeats (e.g., transposable elements or TEs). The presence of tandem repeats was investigated in all contigs using Tandem Repeats Finder (Benson, 1999). TE detection was performed with RepeatMasker² based on the Repbase database³. The presence of snRNA and snoRNA was sought by aligning our genome assembly to the Rfam database⁴ using the cmscan program (Nawrocki and Eddy, 2013). tRNAs were analyzed using tRNAscan-SE (Chan and Lowe, 2019). rRNA was identified as combined Rfam with RNAMMER (Lagesen et al., 2007).

Mycelia of dikaryotic strain 218 were grown on 60-g culture substrate resting in 250-mL glass beakers that had been covered with sterilized polyethylene film and incubated at 32°C. The aerial mycelia that climbed onto the beaker walls (10–14 days) were harvested. The sclerotia were harvested when they reached about 1 cm in diameter.

Total RNA was extracted from each mycelia and sclerotia sample using Trizol reagent, according to the manufacturer’s instructions (TaKaRa, Japan). The quality and quantity of extracted RNA were measured using 1.5% RNAase-free agarose gel electrophoresis, a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, United States), and a Bioanalyzer 2100 system (Agilent Technologies, United States). The construction of cDNA library and RNA-seq of the eligible six RNA samples were carried out by Frasergen, Inc (Wuhan, China) using Illumina HiSeq sequencing.

The sequenced reads were aligned to our sequenced P. tuber-regium genome using Hisat2 (Kim et al., 2015) with default parameters. The aligned records from the aligners in the BAM/SAM format were further checked to remove potential duplicate molecules. Subsequently, gene expression levels were estimated using FPKM values (fragments per kilobase of exon per million fragments mapped) by Cufflinks software (Trapnell et al., 2010).

The EdgeR (v3.6.8) package method was used to evaluate differential gene expression between mycelial and sclerotial samples (Robinson et al., 2010). False discovery rate (FDR) values < 0.05 and | log₂Fold Change| ≥ 1 were used as thresholds to assess the significance of gene expression differences.

Functional analysis of the transcriptome genes was carried out using GO (Ashburner et al., 2000; Young et al., 2010) and KEGG pathways (Kanehisa et al., 2004).

Reverse transcription quantitative polymerase chain reaction (RT-qPCR) was employed to analyze the eight upregulated and eight downregulated candidate genes among transcriptome samples to verify the quality of RNA-seq results. The gene-specific primers used for RT-qPCR are shown in Supplementary Table 1. The reactions were executed in QuantStudio™ 6 Flex System (Thermo Fisher Scientific, United States) using the following parameters: 95°C for 10 min, followed by 40 cycles of 95°C for 15 s, and 60°C for 1 min, and the melt curve for each product was detected. Three biological and three technical replicates were performed for each sample using the combination of E3upl and Tif-5a as the internal reference genes (Sun et al., 2021). The log₂Fold Change of sclerotia relative to mycelia was calculated using the 2^–ΔΔCt method.

In total, the final genome assembly of P. tuber-regium monokaryotic strain 18 yielded a total size of 35.82 Mb, 65 contigs, 47.5% G + C content, and a contig N50 of 4.29 Mb. Approximately, 527 × coverage data were acquired. The monokaryon genome is bigger in size than the dikaryon genome of this previously sequenced mushroom using the Illumina platform (Lam et al., 2018). The genome sequenced in our study is similar in size to the genomes of several other Agaricales species, including L. rhinocerotis (34.31 Mb) (Yap et al., 2014) and P. ostreatus (34.9 Mb) (Qu et al., 2016) but smaller than that of W. cocos (50.48 Mb) (Floudas et al., 2012), P. tuoliensis (48.2 Mb), and P. eryngii (49.9 Mb) (Zhang Z. et al., 2018). The general statistics and view of the sequenced P. tuber-regium genome are shown in Table 1 and Figure 1, respectively.

A large portion of ncRNAs are functional and can produce various regulatory activities, such as signal molecules, ligands, and enzymes, and they are widely expressed in both prokaryotes and eukaryotes (Fu, 2014). The structural ncRNAs, including rRNA, tRNA, snRNA, and snoRNA, and regulatory ncRNAs were identified and mapped to the P. tuber-regium genome (Figure 1). Splicing catalysis by spliceosome is an essential step in the process of gene expression. The spliceosome is formed by the stepwise integration of five snRNPs composed of U1, U2, U4, U5, and U6 snRNAs and more than 150 proteins binding sequentially to pre-mRNA (Stark and Lührmann, 2006). Besides U1, the U2, U4, U5, and U6 snRNAs were also identified in our fungal strain. It is worth noting that U1 was also not identified in other basidiomycetous mushrooms, including P. ostreatus (Qu et al., 2016). With regard to tRNAs, 173 were annotated along with their loci, lengths, and transported amino acids. In addition, six snoRNAs, five 5s rRNAs, five 5.8s rRNAs, four 18s rRNAs, and four 28s rRNAs were identified in the genome.

Repetitive elements are essential for organismal development, including genome organization, evolution, epigenetic diversity, and plasticity (Foulongne-Oriol et al., 2013; Freire-Benéitez et al., 2016). Tandem repeats accounted for 1.13% of the genome assembly with an average length of 127 bp. The same as Agaricus bisporus, microsatellites appeared widely and distributed evenly over the whole genome sequence (Figure 1). The TEs, composed of Class I retrotransposons and Class II DNA transposons, accounted for 2.45% of the genome assembly with a total length of 876,249 bp. For class I, 339 sequences had homologies with long terminal repeats (LTRs), two sequences were observed having long interspersed nuclear elements (LINEs), and none had short interspersed elements (SINEs). For class II, 58 DNA transposons were identified in the genome assembly of P. tuber-regium.

To better elucidate the molecular mechanisms of sclerotium formation in P. tuber-regium, RNA-seq was performed to compare gene expression differences between its mycelia and newly formed sclerotia. After data filtering and trimming, high-quality reads ranging from 29 to 37 Mb (mycelia) and 23 to 31 Mb (sclerotia) were obtained from the six examined libraries (Table 2). The Q30 of the six samples was more than 92%, indicating that the overall sequencing quality was satisfactory, and the data were suitable for further assembly analysis. Hisat2 software mapped over 89.47% of our clean reads to the genome, which was used for further analysis (Table 2). Annotation and functional classification results of 9,815 genes were obtained using the eight data sets. A total of 9,781 (77%) genes associated with known proteins in the Nr database, followed by TrEMBL database (9,769: 77%), Pfam database (6,623: 52%), SwissProt database (5,318: 41%), KOG database (4,264: 33%), GO database (4,210: 33%), KEGG database (2,922: 23%), and COG database (2,746: 21%).

According to DEG analysis, we found conspicuous differences in the gene expression levels between mycelia and newly formed sclerotia. We compared libraries for gene expression differences that were ≥ 2 or ≤ −2 and detected 1,146 upregulated genes and 1,249 downregulated genes in the newly formed sclerotia compared with mycelia. To verify the reliability of the transcriptomic data, RT-qPCR of 16 DEGs that had been observed through RNA-seq showed that our profiling of these DEGs were consistent with RNA-seq results (Figure 2).

GO enrichment was conducted for the genes obtained from each group of samples, which were classified into the three primary GO categories annotated to biological process (GO-BP), cellular component (GO-CC), and molecular function (GO-MF) (Figure 3A). As shown from the results of GO-BP analysis, the most highly enriched term was metabolic processes, followed by cellular processes, illustrating that these processes are of functional importance for mycelial growth and sclerotial development in P. tuber-regium. In the category of GO-CC, membrane and cell were significantly enriched in the transcriptome genes of mycelia and sclerotia. In the GO-MF category, the most highly enriched term was binding. The second highly enriched GO-MF category was catalytic activity. Overall, the GO terms for the transcriptome genes of mycelia and sclerotia, which included metabolic process, membrane, and catalytic activity might potentially account for the formation of sclerotia in P. tuber-regium.

For better understanding of the functional pathways involved in sclerotia differentiation, we matched the transcriptome genes of mycelia and sclerotia to the KEGG database and displayed several metabolic and signal pathways (Figure 3B). We also conducted the KEGG enrichment analysis of DEGs, the KEGG results from our two groups of libraries revealed that DEGs listed in differentiation of mycelia to sclerotia were primarily enriched in the carbon metabolism; lysosome, cysteine, and methionine metabolism; and arginine and proline metabolism, etc (p < 0.05).

Some studies revealed that a series of metabolic process participated in sclerotium formation (Yap et al., 2015; Zhang et al., 2016). The analysis of these DEGs mainly focused on carbohydrate metabolism, oxidation-reduction process, and cell wall synthesis. In addition, some gene clusters speculated to participate in the sclerotia development were further analyzed, such as carbohydrate active enzymes (CAZymes), oxidoreductases, chitin synthases, and hydrophobins.

The differential expression of several genes in the CAZymes correlated well with our transcriptomic data. The CAZymes, including glycoside hydrolases (GHs), carbohydrate esterases (CEs), auxiliary activities (AAs), carbohydrate-binding modules (CBMs), glycosyl transferases (GTs), and polysaccharide lyases (PLs), are involved in carbohydrate metabolism. We analyzed the DEGs between mycelia and sclerotia and identified 76 and 97 CAZymes genes that were significantly upregulated and downregulated during the early stage of sclerotial development, respectively (Figure 4). Among these DEGs, more than half of them belonged to GH (33.53%) and AA (30.64%) families, indicating the importance of P. tuber-regium GHs and AAs families for sclerotium formation. In CAZymes, the CBMs encoded by maker-contig_6_pilon-exonerate_protein2genome-gene-13.1 and evm.TU.contig_17_pilon.35 were the most upregulated and downregulated in scletotia, respectively. It was reported that the CAZymes genes were not only involved in substrate decomposition, but also in sclerotium formation in W. cocos (Zhang et al., 2016). Our findings suggested that the CAZymes genes identified in P. tuber-regium might also participate in sclerotium formation.

We analyzed the D-glucose, the monomer of glucan and 1, 3-ß-glucan related pathway derived from the starch and sucrose metabolism pathway in KEGG. The beta-glucosidase genes (evm.TU.contig_30_pilon.59 and evm.TU.contig_5_pilon.88) and glucan 1, 3-beta-glucosidase genes (evm.TU.contig_3_pilon.1871 and maker-contig_32_pilon-exonerate_protein2genome-gene-7.3) involved in the D-glucose pathway were downregulated with the exception of gene maker-contig_6_pilon-exonerate_protein2genome-gene-7.10. This may account for the fact that sclerotia may only have early stage growth needs and typically exist in relatively dormant states during other growth periods when requirements of energy for life cycle activities are reduced. Additionally, with the exception of evm.TU.contig_93_pilon.104, most of the chitinase genes (evm.TU.contig_13_pilon.236, evm.TU.contig_15_pilon.689, and evm.TU.contig_3_pilon.1375) were downregulated in our P. tuber-regium sclerotia.

Genes related to the oxidation-reduction process were differentially expressed during sclerotium formation in P. tuber-regium. Oxidation -reduction–related genes were differentially expressed in mycelium vs. sclerotium production in other mushroom-producing species, such as Polyporus umbellatus, W. cocos, and L. rhinocerotis (Song et al., 2014; Yap et al., 2015; Wu et al., 2016; Li et al., 2017). The molecular mechanisms controlling sclerotium formation might include oxidative stress, signal transduction, and gene expression regulation (Sun et al., 2020). In this study, 155 genes involved in oxidation-reduction progress were significantly and differentially expressed in GO-BP (Figure 5A). Among them, 31 Cytochrome P450 (CYP) monooxygenases, demonstrating differential transcriptional abundances, deserved further attention. CYP monooxygenases have the unique ability to catalyze regio-, chemo-, and stereospecific oxidation of a wide range of substrates under mild reaction conditions, thereby addressing a significant challenge in chemocatalysis (Durairaj et al., 2016). Among the DEGs, the genes (evm.TU.contig_17_pilon.68 and evm.TU.contig_17_pilon.69) encoding laccases were dramatically upregulated in sclerotia. Compared with the mycelia stage, the lipoxygenase encoded by evm.TU.contig_25_pilon.317 in the sclerotia decreased significantly.

Differentially expressed genes involved in the cell wall system might participate in sclerotium formation in P. tuber-regium. Cell wall biogenesis is a very complicated process that incorporates not only multiple gene families directly involved in the polysaccharide biosynthesis and the remodeling of different cell wall components, but also structural proteins involved in cell wall assembly and rearrangement (Hu et al., 2017). We further analyzed genes with upregulated or downregulated expression in the formation of sclerotia that associated with GO cell composition terms related to “wall.” Notably, 16 cell wall–associated genes were screened, and 14 of them were hydrophobins (HPs), the structural constituents of cell walls. HPs are surface-active proteins unique in fungal cells and have roles in fungal growth as structural components and in the interaction of fungi with their environment (Linder et al., 2005; Kulkarni et al., 2017). In P. tuber-regium sclerotia, evm.TU.contig_38_pilon.140 had the highest expression abundance, with more than 50,000 FPKM (Figure 5B). In particular, expression level of HP gene evm.TU.contig_6_pilon.674 in sclerotia was hundreds of times greater than in mycelia. It is reported that HPs can affect the cell wall composition and protect emergent structures against adverse environmental conditions (Wetter et al., 2000; Wösten, 2001). We also observed an upregulated PII-type proteinase homologous to Lactococcus lactis subsp. Cremoris and a downregulated minor extracellular protease vpr homologous to Bacillus subtilis (strain 168) related to serine-type endopeptidase activity. Cell walls of sclerotia were thicker than those of mycelia in P. tuber-regium (Sun et al., 2020). This is likely because of the accumulation of cell wall components in addition to HPs, such as chitin, glycan, and xylan, as reported for P. umbellatus sclerotia (Li et al., 2021).