The mycelium of C. lacerata CGMCC No. 10485 was stored at 4°C after initial cultivation on PDA solid medium for 7 days at 25°C. A small amount of the mycelium from slant culture was transferred into a 500-mL conical shake flask containing 150 mL first-order seed culture medium (soluble starch 20 g/L, spray-dried corn steep liquor 6 g/L, KH₂PO₄ 1 g/L, sterilized at 121°C for 20 min), then cultivated for 3.5 days in a rotary shaker at 150 rpm and 25°C.

Then, 7.5 mL of the first-order seed culture was used to inoculate a 500-mL flask containing 150 mL of second-order seed culture medium (glucose 80 g/L, spray-dried corn steep liquor 8 g/L, KH₂PO₄ 5 g/L, sterilized at 115°C for 30 min) and grown in a rotary shaker at 150 rpm and 25°C for sequential sampling on the 2nd day (lag phase, denoted as CL1), day 2 (early logarithmic phase, CL2), day 5 day (logarithmic phase, CL3), and day 8 (stationary phase, CL4) (Growth characteristics were described in Supplementary Figure 1). The mycelium pellets from different culture stages were separated from the medium by filtration, shock-frozen in liquid nitrogen, and stored at −80°C for RNA/DNA isolation and sequencing.

The genomic DNA was isolated using the DNeasy Plant Mini Kit (QIAGEN) according to the manufacturer’s instruction. PacBio sequencing was performed by Novogene (Tianjin, China) on a PacBio RS platform with one SMART cell, which generated long-read data. In addition, another run of Illumina sequencing was performed on a NovaSeq 6000 platform (standard 2 × 150 paired-end libraries), which generated short-read data. All the raw data in FASTQ format were processed to get clean reads by removing adapters, N bases, and low-quality bases.

Total RNA collected at CL1-CL2-CL3-CL4 was extracted using the RNAprep Pure Cell/Bacteria Kit (TIANGEN, China). The integrity of the total RNA was assessed using 1.0% agarose gel electrophoresis. The concentration and purity were determined using a NanoDrop 2000 spectrophotometer. Total RNA (1 μg) was treated with the QIAseq™ FastSelect−5S/16S/23S Kit (QIAGEN) to remove rRNA, and cDNA libraries were prepared using the QIAseq™ Stranded Total RNA Lib Kit according to the manufacturer’s instructions. Illumina sequencing was performed on a NovaSeq 6000 platform by Novogene (Tianjin, China). Raw data were filtered to obtain the high-quality RNA-seq data using the same procedure described above for genomic short-read data. The raw sequencing data (both genomic and transcriptomic) have been submitted to the NCBI Sequence Read Archive under BioProject PRJNA804482 (link).¹

Long-read data generated by the PacBio platform were de novo assembled using CANU 2.1 as described previously (Koren et al., 2017). Briefly, the long reads were corrected, trimmed, and then assembled into contigs. Subsequently, the contigs were polished based on Illumina short reads using Pilon 1.24 (Walker et al., 2014). The genome completeness was evaluated using BUSCO (version 3.0.2) with fungi_odb10 OrthoDB database (Simao et al., 2015).

The RepeatMasker program (version 4.1.1)² was used for masking repetitive elements in nucleotide sequences of contigs. The genomic structure annotation was performed using BRAKER2 (version 2.1.6) (Bruna et al., 2021), which allows for fully automated training of the gene prediction tools GeneMark-ET (Lomsadze et al., 2014) and AUGUSTUS (Stanke et al., 2006) using RNA-seq data. Microsatellite identification tool (MISA) (Beier et al., 2017) was used to determine the distribution and frequency of various types of simple sequence repeats (SSRs). tRNAs and rRNAs were identified using tRNAscan-SE (version 2.0.7) (Chan and Lowe, 2019) and RNAmmer (version 1.2) (Lagesen et al., 2007), respectively. Other non-coding RNAs, including small RNAs (sRNAs) and small nuclear RNA (snRNAs), were inferred from Rfam using Infernal (version 1.1.4) (Nawrocki and Eddy, 2013).

Functional annotation of genes was conducted by homology searching against protein sequences from the SwissProt,³ NR⁴ and eukaryotic orthologous groups (KOG)⁵ databases using Diamond 2.0.9 (Buchfink et al., 2015). For multiple matches of a single protein, only the best match was retained. The gene ontology (GO) annotation was carried out using InterProScan (version 5.45) (Jones et al., 2014). The functional pathway annotation (Kyoto Encyclopedia of Genes and Genomes, KEGG) was carried out using BlastKOALA in KEGG (Kanehisa et al., 2016). Signal peptides and cleavage sites of C. lacerata proteins were predicted using SignalP (version 4.1) (Almagro Armenteros et al., 2019). All proteins with signal peptides were analyzed for the presence of transmembrane domains using TMHMM (version 2.0) (Moller et al., 2001). Similarly, Phobius (version 1.01) (Kall et al., 2004),⁶ which predicted the transmembrane structure and signal peptide, was used to expand the subsequent secretory proteins. A protein containing a signal peptide and no transmembrane domain was identified as a secretory protein. The genes encoding Carbohydrate-Active enZYme (CAZyme) (Lombard et al., 2014) were annotated using dbCAN2 (Zhang et al., 2018). The protein sequences of map00900, map00902, map00904, and map00909, which are the major biosynthetic pathways involved in terpenoid synthesis, were downloaded from the KEGG database to evaluate the presence/absence of terpenoid biosynthesis genes.

The 18 representative genomes of species from the order Polyporales containing both sequences and genomic annotation information were downloaded from NCBI and listed in Table 1. Orthofinder (Emms and Kelly, 2019) and FigTree v1.4.4⁷ were used to compute orthogroups and construct genome-level phylogenetic trees with single-copy orthologous genes based on Maximum Likelihood (ML) algorithm. FastANI (Jain et al., 2018) was used to calculate the whole-genome average nucleotide identity (ANI) with the parameter “fragLen” being set to 1000. Additionally, ITS sequences (Supplementary Material) from 61 representative TYPE materials (Federhen, 2015) of these 18 species were collected from the NCBI fungus ITS project for phylogenetic analysis in MEGA 7 (Kumar et al., 2016), using Neighbor-Joining method with default parameters.

The C. lacerata genome sequencing yielded 16 Gb of PacBio RS long-read data (∼400 × coverage), 20 Gb of Illumina short-read data (∼500 × coverage), and 12 Gb of Illumina cDNA data for the hybrid assembly and genome annotation. The genome was initially assembled using Canu via reads correction, trimming and assembly based on PacBio long-fragment sequences, then polished using Pilon based on Illumina short reads. The final C. lacerata genome assembly contained 58 scaffolds with a total consensus genome size of 36,361,585 bp (∼36 Mb) and GC content of 49.33%, which was comparable to the genome sizes of previously sequenced fungi from the order Polyporales (28∼60 Mb). The maximal scaffold size was 4,415,373 bp, and the N50 scaffold size was 3,409.20 kb. The BUSCO evaluation showed 98.4% completeness of this genome, indicating a high-quality genome assembly (Table 2). All annotation statistics are listed in Supplementary Table 1.

A total of 13,243 genes were predicted in the C. lacerata genome with an average length of 1,860 bp, 9,085 of which were successfully annotated using the NR, SwissProt, and KOG databases. In addition, 223 non-coding RNAs, including tRNAs, rRNAs, and snRNA were identified in the C. lacerata genome. Other information, such as the numbers of exons and introns, is listed in Table 2.

Non-coding repeat regions such as microsatellites and simple sequence repeats (SSRs) provide excellent information for assessing genomic variation and are thus frequently used as molecular markers for distinguishing even closely related strains. A total of 5,475 repetitive sequences (277,515 bp) were identified, accounting for 0.76% of the whole genome. More than 81% of these were simple sequence repeats, followed by transposable elements (TEs), including LTRs, LINEs, SINEs, and transposons (Figure 1A). Moreover, 2,061 SSRs were found in the genome (Figure 1B), nearly 54% of which were mononucleotide tandem repeats, followed by tri-, di-, tetra, penta-, and hexanucleotide motifs.

Overall, 6,083 genes (45.9% of the total protein-coding genes) were annotated with GO terms classified as biological process (3,828 genes), molecular function (5,346 genes), and cellular component (1,544 genes), respectively. The top GO terms were protein binding, oxidation-reduction process, ATP binding, nucleic acid binding, zinc ion binding, integral component of the membrane, protein phosphorylation, membrane, regulation of transcription, DNA-templated, carbohydrate metabolic process, nucleus, and host cell nucleus (Figure 2A). In the KEGG annotation of C. lacerata genome, 4,215 genes were successfully annotated with 3,170 KEGG Orthologous (KO) terms. These terms were grouped into 42 KEGG BRITE functional hierarchies, including enzymes, membrane trafficking, chromosome and associated proteins, ribosome biogenesis, mitochondrial biogenesis, messenger RNA biogenesis, exosome, DNA repair and recombination proteins, ubiquitin system, and spliceosome (Figure 2B). Both GO and KEGG annotations were mainly concentrated in the basic metabolism and host interactions on which C. lacerata depends.

The KOG database is a eukaryote-specific version of the COG (clusters of orthologous groups of proteins) database. Proteins orthologs typically occupy the same functional niche in different species and refer to proteins evolved from vertical families (speciation) from different species. In this study, a total of 2,721 genes (19.0% of the total) were annotated with KOG terms and grouped into 23 classes, generally covering the essential functions of metabolism, genetic information processing, environmental responses, and cellular processes. The top KOG terms were R: General function prediction only, O: Posttranslational modification, protein turnover, chaperones, J: Translation, ribosomal structure and biogenesis, C: Energy production and conversion, as well as A: RNA processing and modification (Figure 2C).

The “secretome” constitutes the entire set of secreted proteins of a microorganism (Ganesan, 2016), including functionally diverse classes of molecules, such as digestive enzymes, chemokines, antibodies, extracellular proteinases, morphogens, toxins, and antimicrobial peptides. For C. lacerata, 1,163 genes were found to encode secretory proteins that contain N-terminal signal peptides and no transmembrane domains, accounting for 8.8% of all proteins across the whole genome. Among them, 527 were annotated with GO terms (Figure 2D). Significantly, the top-ranked terms covered a number of secretory digestive enzymes required by this fungus to colonize the woody host and conduct its saprophytic lifestyle, such as protein binding (GO:0005515, protein degradation tagging, cell adhesion), hydrolyzing O-glycosyl compounds (GO:0004553, catalysis of the hydrolysis of any O-glycosyl bonds for cellulose degradation), aspartic-type endopeptidase activity (GO:0004190), serine-type carboxypeptidase activity (GO:0004185), and proteolysis, serine-type endopeptidase activity (GO:0006508), and metalloendopeptidase activity (GO:0004222). Additionally, a few secretory proteins were involved in a host of diverse and vital biological processes, including calcium ion binding (GO:0005509), zinc ion binding (GO:0008270), signal transduction (GO:0007165), cell surface receptor signaling pathway (GO:0007166), mycotoxin biosynthetic process (GO:0043386), and regulation of cell cycle (GO:0051726), which are related to signal transduction, proliferation, survival, defense, and virulence factors.

Even though several studies have been conducted on the phylogeny of C. lacerata, most of them were performed using ITS or 18S rDNA sequence alignment (Jia et al., 2013; Miettinen et al., 2016; Ponlada et al., 2016; Wu et al., 2017; Yuan et al., 2017; Chen et al., 2020). In the present study, to reveal the primary phylogenetic position of C. lacerata CGMCC No. 10485, a phylogenetic tree was reconstructed (Figure 3A) with 62 ITS sequences (Federhen, 2015) associated with the 19 strains. Collectively, the ITS samples from Ceriporia, Phlebia, and Phanerochaete were more closely clustered than any other species from the order Polyporales, which was similar to previous reports (Wu et al., 2017; Yuan et al., 2017; Chen et al., 2020). Evidently, Ceriporia strains formed a clear monophyletic group, indicating the taxonomic and phylogenetic position of C. lacerata CGMCC No. 10485. Subsequently, a denser sampling of larger and orthologous gene sets across the whole genome was conducted to construct more detailed phylogenomic trees to advance the phylogeny and systematics of the order Polyporales.

Genome-based phylogenetic analyses, such as multiple alignment of single-copy orthologous genes (SCOG) and average nucleotide identity (ANI), provide a more robust and reliable taxonomy and phylogeny due to their greater resolution of evolutionary relatedness between different species (Zhang et al., 2017). As shown in Figure 3B, the SCOG phylogenetic tree showed that Phlebia, Ceriporia, Phlebiopsis, and Phanerochaete were clustered and merged into an evolutionary branch. ANI is a measure of nucleotide-level genomic similarity, for both complete and draft assemblies, and is a superior indicator of relatedness between the coding regions of two genomes. Variable evolutionary rates as well as the horizontal gene transfer of one or a few genes have no influence on ANI values because the effect of fast-evolving genes is balanced by the slow evolution of other genes (Jain et al., 2018). Generally, an ANI value greater than 95% is usually used as the gold standard for species classification and clustering. According to the results shown in Figure 3C, the ANI values variated between 70 and 80%. This was significantly lower than 95%, confirming that these genomes were from different species. ANI-based clustering indicated that the phylogenetic relatedness between any two genomes was consistent with the SCOG analysis, which further emphasized the robustness and reliability of the phylogenomic analysis shown in Figure 3B.

Carbohydrate-active enzymes (CAZymes) that can act synergistically on a wide range of glycosidic monomers, oligomers, or polymers in a regio- or stereo-specific manner, are required for the breakdown of carbohydrates, along with lignin and hemicellulose. In general, CAZymes include glycoside hydrolases (GHs), glycosyltransferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), auxiliary activity proteins (AAs), and carbohydrate-binding modules (CBMs).⁸

A total of 806 CAZyme-encoding genes were identified in the genome of C. lacerata, encompassing 41.8% GHs, 21.8% GTs, 1.74% PLs, 7.32% CEs, 13.5% AAs, and 13.8% CBMs (Figure 4A). Notably, phylogenetic analysis based on CAZyme sequences (the left pane of Figure 4A) was in excellent agreement with the results of SCOG and ANI in Figures 3B,C, indicating the species-specific evolution of CAZymes among Polyporales for their unique saprophytic lifestyle growing on a wood-based matrix.

To evaluate the CAZYme of C. lacerata CGMCC No. 10485, comparisons with other 18 fungi of Polyporales were performed. As shown in Figure 4B, CAZymes can be further categorized into subfamilies based on amino acid sequence similarity, substrate selectivity, and catalytic mechanism (Cantarel et al., 2009). Endoglucanases (EC 3.2.1.4), exoglucanases (syn: cellobiohydrolases; EC 3.2.1.91), and β-glucosidases (EC 3.2.1.21) produced by white-rot fungi, are able to break down cellulose synergistically (Daniel, 2016). The endoglucanases in the CAZyme database mainly belong to the GH5∼10 families, while others are classified into GH12, 26, 44, 45, 48, 51, 74, 124, and 148. In the genome of C. lacerata, several genes that encode enzymes of GH5, GH7, GH9, GH10, GH12 and GH45 were annotated, among which GH5 (including GH5_5, GH5_50 and GH5_7) and GH9 were more abundant compared with the other species used in this study (Figure 4B). Apart from endoglucanases, several exoglucanases, including GH5, 9, and 51, and β-glucosidases, GH1, 2, 3, 5, 16, and 30 were successfully annotated, demonstrating the potential of this strain to degrade cellulose.

Compared to a wide variety of cellulases with different functions, white rot ligninolytic enzymes are mainly grouped into two families: AA1 [laccase/p-diphenol: oxygen oxidoreductase/ferroxidase (EC 1.10.3.2)] and AA2 [manganese peroxidase (EC 1.11.1.13), lignin peroxidase (EC 1.11.1.14), versatile peroxidases (EC 1.11.1.16)] (Mkel et al., 2021). Lignin peroxidases (EC 1.11.1.14) are powerful oxidants with a high redox potential that oxidize lignin’s non-aromatic components. Manganese peroxidase (EC 1.11.1.13) is a Mn-dependent enzyme that can oxidize aromatic substrates but is inactive on non-aromatic lignin. The catalytic activity of manganese and lignin peroxidases are combined in versatile enzymes (EC 1.11.1.16). Laccases (EC 1.10.3.2) are members of the multi-copper oxidase family that catalyze a one-electron oxidation coupled with a four-electron reduction of molecular oxygen to water (Abdel-Hamid et al., 2013). In the genome of C. lacerata, 5 and 12 genes were identified as encoding enzymes of AA1 and AA2, respectively. We also discovered a large number of additional AA genes in the genome of C. lacerata. They were mainly lytic polysaccharide monooxygenases and enzymes involved in the biodegradation of aromatic compounds, including AA9, AA3_ 2, AA5_ 1, AA10, AA6, AA3_ 3, AA7, AA14, and AA4. Among them, AA4 and AA9 had the greatest gene counts among the selected species (Figure 4B). AA4 proteins are vanillyl-alcohol oxidases that catalyze the conversion of a wide range of aromatic compounds bearing side chains at the para-position of the aromatic ring, while AA9 proteins are copper-dependent lytic polysaccharide monooxygenases responsible for the cleavage of cellulose chains through oxidation of carbons C1, C6, and/or C4. The diversity and abundance of AA families in the C. lacerata genome reflects its strong capacity to degrade lignin and various PACs (Mori and Kondo, 2002; Suhara et al., 2003; Xu et al., 2008; Lin et al., 2011; Hong et al., 2012, 2013; Wang et al., 2017; Yanto et al., 2019; Al-Hawash et al., 2020).

Terpenoids, also known as isoprenoids, are categorized by the number and arrangement of carbon atoms generated by the linear arrangement of isoprene units followed by cyclization and rearrangements of the carbon skeleton. The terpenoid backbone is biosynthesized via the mevalonic acid (MVA) pathway and the 2-C-methyl-d-erythritol 4-phosphate/1-deoxy-d-xylulose 5-phosphate (MEP/DOXP) pathway (Figure 5A). In the KEGG database, map00900, map00902, map00904 and map00909 are the major biosynthetic pathways that synthesize the backbone of monoterpenoids, diterpenoids, sesquiterpenoids and triterpenoids (Figure 5A). To evaluate the terpenoid synthesis potential of C. lacerata CGMCC No. 10485, another 8 genomes, including Phlebiopsis gigantea, Phlebia centrifuga and Phanerochaete carnosa (with higher relatedness), as well as the typical terpenoid-producing species Ganoderma sinense, Trametes versicolor, T. pubescens, T. coccinea, and T. cinnabarina (Ludwiczuk et al., 2017; Winska et al., 2019), were analyzed to annotate their terpenoid biosynthetic pathways. The EC numbers of each gene mapped onto the pathways were converted into KO identifiers.

The annotated KO terms from each pathway are summarized in Figures 5B–E. A total of 131 KO terms involved in terpenoid synthesis were found in the genome of C. lacerata. It’s worth noting that all the KOs in each pathway could be classified into two groups based on their presence or absence in the genomes, characteristic KOs that were shared among nine species, and fingerprint KOs that were unique to one or several species.

The MVA and MEP/DOXP pathways make up the main part of map00900. Among the nine genomes, 53∼59 KO terms can be found (Supplementary Table 3). Interestingly, the characteristic KOs were primarily found in the MVA pathway (K00626, K01641, K00021, K00869, K00938, and Figure 5B), farnesyl diphosphate synthesis pathway (K15890, K15891, K05906, K15889, K00587, K06013, and K08658), and geranylgeranyl diphosphate pathway (K05355). However, the absence of downstream strain-specific KOs (K00099, K12506, or K00919) resulted in an incomplete MEP/DOXP pathway, forming the fingerprint area to distinguish the nine species. This finding demonstrated that the MVA pathway was the main determinant of terpenoid synthesis among these white-rot fungi.

According to the KO annotations shown in Figure 5C (map00902, Supplementary Table 4), 22 ∼ 27 KO terms can be found for monoterpenoid biosynthesis pathway. Unique to C. lacerata were K21069, K07381, K05525, and K15090, which are predicted to be involved in the biosynthesis of (R)-ipsdienol, (6E)-8-hydroxylinalool, trans-8-oxolinalool, and (-)-trans-carveol, respectively.⁹ In pathway 00904 (Figure 5D and Supplementary Table 5), 20∼ 26 KO terms were annotated. The fingerprint KOs of C. lacerata were mainly enriched in the synthetic pathways of ferruginol, 11-hydroxyferruginol and salviol (diterpenoids).¹⁰ In pathway 00909 (Figure 5E and Supplementary Table 6), 22 ∼ 27 KOs were identified in the nine genomes. Characteristic KOs, K00801 and K00501 are the critical enzymes that successively convert farnesyl diphosphate into triterpenoids. Sesquiterpenes are produced from farnesyl diphosphate via one- or two-step catalytic reactions catalyzed by strain-specific enzymes. In this pathway, fingerprint KOs encode various enzymes that depend on the species. In C. lacerata, the fingerprint KOs were mainly related to the biosynthesis of nerolidol (sesquiterpene) and germacrene A (triterpenoid).¹¹

In conclusion, the characteristic and fingerprint KOs provide an effective identifier to distinguish the main or branching KEGG pathways and are of great importance for understanding strain-specific terpenoid synthesis enzymes.