The monokaryotic strain was isolated from a wild fruiting body collected from Changbai Mt., Antu County, Yanbian, Jilin Province, China. The mycelia of D. sinensis were harvested after growing on sterile cellophane covering potato dextrose agar (PDA) medium at 28°C for 10–14 days. The mycelia were collected, snap frozen in liquid nitrogen, and then stored at −80°C for further use. The enriched genomic DNA was extracted from the mycelia with a QIAGEN Genomic Kit and detected by agarose gel electrophoresis.

PacBio Sequel and MGISEQ2000 platforms were used for whole genome sequencing at Wuhan Nextomics Biosciences Co., Ltd. (Wuhan, China).

The quality of the DNA was examined using four methods: the DNA was inspected whether the appearance and shape of the sample contained foreign bodies; degradation and DNA fragment size of samples were detected by 0.75% agarose electrophoresis; DNA purity was measured by NanoDrop spectrophotometer; and precise quantification of the DNA was performed using a Qubit system.

After the samples passed the quality inspection, the genomic DNA was sheared with g-TUBEs (Covaris, United States) according to the size of the fragments built in the library, and magnetic beads were used to enrich and purify the targeted fragments of DNA. The fragmented DNA was then repaired for damage and terminal repair. Stem loop sequencing splices were connected at both ends of the DNA fragment, and exonuclease was used to remove the fragments that failed to connect. Next, BluePippin (Sage Science, United States) was used to screen the target fragments, and the library was purified. The library fragment size was then measured with an Agilent 2,100 Bioanalyzer (Agilent Technologies, United States). After construction of the library, DNA templates and enzyme complexes with a certain concentration and volume were transferred to the nanopore of the PacBio Sequel series sequencer for real-time single molecule sequencing.

Quality control was conducted on the raw sequencing data, and the low-quality areas and adapter sequences were removed to obtain high-quality DNA sequence data. Hifiasm software (parameter: −n5) was used for pure three-generation assembly to obtain the preliminary assembled genome sequence. Nextpolish software was used to calibrate the genome sequence four times with the second-generation data, and the final genome sequence was obtained.

For genome repeat sequence prediction, GMATA version 2.2 (Wang and Wang, 2016) was used to analyze the simple sequence repeats in the genome. TRF version 4.07b (Benson, 1999) was used to analyze tandem repeats in the genome with default parameters. The interspersed repetitive sequences were predicted using RepeatMasker version open-4.0.9 (Bedell et al., 2000). Three methods were employed in gene structure prediction, namely, transcriptome prediction, homologous protein prediction and de novo prediction. Based on transcriptome sequence data, genome gene prediction was performed by PASA version 2.3.3 (Haas et al., 2003). Homologous protein prediction used GeMoMa version 1.6.1 (Keilwagen et al., 2016) to compare the corresponding protein information with the genome. In addition, de novo prediction was performed using Augustus version 3.3.1 (Stanke et al., 2008). Subsequently, the gene prediction results obtained by the above three methods were integrated by using EVidenceModeler version 1.1.1 (Haas et al., 2008). Genome integrity was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO) version 3.0.1 (Simão et al., 2015).

Gene functions were predicted with references to nine databases: Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/) to understand the higher functions and utility of cells, organisms, and ecosystems at the level of molecular information and to reveal the hidden features in biological data (Kanehisa et al., 2017), Gene Ontology (GO) database (http://www.geneontology.org) to understand the living organisms in three pathways: molecular function, cellular component, and biological process (Ashburner et al., 2000); Eukaryotic Orthologous Group of Protein (KOG) database (https://www.creative-proteomics.com/services/kog-annotation-analysis-service.htm), a eukaryote-specific version of the Clusters of Orthologous Groups tool for identifying the ortholog and paralog proteins (Galperin et al., 2015); Non-Redundant Protein (NR) database (https://www.ncbi.nlm.nih.gov/protein/), a non-redundant protein library for protein function annotation; Fungal Cytochrome P450 (CYP) database (http://p450.riceblast.snu.ac.kr/cyp.php) and CYP Engineering database (http://www.cyped.uni-stuttgart.de) to classify and analyze the CYP monooxygenase family for a better understanding of their biochemical characteristics and sequence-structure-function relationships (Sirim et al., 2009; Moktali et al., 2012); Swiss-Prot database (https://www.uniprot.org/), a protein sequence library designed to provide the high levels of annotation, minimal redundancy, and integration with other databases (Stanke and Waack, 2004; Stanke et al., 2006); Pfam database (http://pfam.xfam.org/), an aggregation of protein families including comparative sequences, species information, and hidden Markov models corresponding to special gene families (Mistry et al., 2021); Carbohydrate-Active Enzyme (CAZyme) database (http://www.cazy.org/) (Drula et al., 2022), a special database dedicated to presentation and analysis of genomic, structural, and biochemical information on CAZymes. CAZyme database classifies CAZymes into different protein families based on the similarity of amino acid sequences in the protein domain, which covers related enzymes required for lignocellulosic degradation. At present, the CAZyme database is an important reference for annotating CAZymes in most microbial genomes or metagenomes (Cantarel et al., 2009; Zhu et al., 2016; Busk et al., 2017; Helbert et al., 2019).

All predicted coding genes were aligned with these nine databases with E-value cut-offs ≤1 × 10⁻⁵, identity ≥40%, and coverage >30%.

The gene clusters of secondary metabolites were predicted using antiSMASH version 7.0.1 with default parameters (Blin et al., 2023).

To explore the evolutionary dynamics of D. sinensis, the genome sequences of 16 additional fungal species were downloaded from the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/genbank/) for phylogenomic analysis.

Single-copy orthologous genes from the 17 fungal species were inferred using OrthoFinder version 2.5.4 (Emms and Kelly, 2019) with the mafft option for subsequent multiple sequence alignment. On the basis of the resulting alignment, a maximum-likelihood tree was reconstructed using RAxML version 8.2.12 (Stamatakis et al., 2008). In addition, Timetree was used to calculate the divergence time among the 17 fungal species (Kumar et al., 2022).

Expansion and contraction of gene families were determined using CAFE version 4.2.1 (Han et al., 2013) with the following parameters: a cut-off p-value of 0.05; number of random samples = 1,000; the lambda value to calculate birth and death rates.

The internal transcribed spacer (ITS) sequences in Daedaleopsis species were aligned using ClustalX 1.83 (Chenna et al., 2003) and optimized manually using BioEdit 7.0.5.3 (Hall, 1999) prior to phylogenetic analysis. Thereafter, maximum parsimony (MP) bootstrap analysis was performed using PAUP^* version 4.0 beta 10 to reveal the phylogenetics of Daedaleopsis species (Swofford, 2002). ITS sequence of D. sinensis was extracted as 559 bp from its genome data and then compared with other ten species in the same genus. At the same time, Earliella scabrosa was selected as an outgroup to construct the evolutionary tree (Li et al., 2016).

To explore the dynamics of speciation in D. sinensis, the genome sequences of D. sinensis and D. nitida were aligned in pairs using MCScanX (Wang et al., 2012) and NGenomeSyn (He et al., 2023).

In addition, the numbers of genes encoding various families of CAZymes in D. sinensis and D. nitida were clustered in a heatmap using TBtools version 1.127 (Chen et al., 2020) with the log-scale option.

Strain was identified as D. sinensis based on internal transcribed spacer-barcoding sequences. Figure 1 shows that the species used in the present study clustered with D. sinensis JX569732 and KU892444. Only one base difference from D. sinensis JX569732 and two base differences from D. sinensis KU892444, further confirming the strain used is D. sinensis based on phylogenetics.

The 51.67 Mb genome sequence was assembled from 6,113 Mb raw data (106 × genome coverage) and to 24 contigs with a GC content of 56.46% (Table 1). Of the 24 contigs, the longest one was 5.41 Mb, whereas the N50 length was 4.40 Mb (Table 2). The GC skew did not present an obvious distribution pattern in the whole genome (Figure 2). A K-mer analysis with a depth of 34 showed a 1.60% heterozygous rate of the genome. Collectively, the above information indicates the high quality of the genome sequence assembly of D. sinensis.

A total of 12,153 protein-coding genes were predicted with an average gene length of 2,420.31 bp and an average coding sequences (CDS) length of 1,476.57 bp, whereas 280 non-coding RNAs (ncRNAs) were predicted, accounting for 23.93% of the whole genome sequence. The total length of repetitive sequences was 7.62 Mb (Table 2) with the tandem repeat sequences and the interspersed repeated sequences accounting for 7.73% and 12.11%, respectively, of the whole genome sequence. The length of the transposable elements was 4.60 Mb with long terminal repeats (LTR) and non-LTR accounting for 7.71% and 1.19%, respectively, of the whole genome sequence. Assessment of genome integrity using BUSCO showed 97.36% of complete genome, indicating that the vast majority of conserved genes were predicted to be relatively complete and implying that there was high confidence in the genome assembly and prediction.

The annotation of genes using the KEGG pathway database can expand understanding of the biological functions and interactions of genes. A total of 3,831 annotated genes had matched in the KEGG database and were assigned into five levels. Of the five main levels, metabolism was the largest (2,400 genes, 48.40%), followed by genetic information processing (826 genes, 16.66%), organismal systems (737 genes, 14.86%), cellular processes (665 genes, 13.41%), and environmental information processing (331 genes, 6.67%). Thus, there were genes corresponding to significant metabolic processes. Many pathways related to lignocellulosic degradation and some key enzymes secreted by D. sinensis were identified in the KEGG database.

Regarding the metabolism and biosynthesis of sugar-containing compounds, 22 enzymes encoded by 43 genes and involved in the biosynthesis of polysaccharides (starch and sucrose metabolism) were identified from D. sinensis (Supplementary Figure S1; Supplementary Table S1). Most of these enzymes were encoded by single- or double-copy genes, whereas β-glucosidase and cellulose 1,4-β-cellobiosidase were encoded by eight- and four-copy genes, respectively. In the glycolysis/gluconeogenesis pathway (Supplementary Figure S2; Supplementary Table S2), there were 24 enzymes encoded by 35 genes, some of which had double or multiple copies. Seventeen key enzymes involved in terpenoid backbone biosynthesis via the mevalonate pathway were identified from D. sinensis (Supplementary Figure S3). Of these 17 key enzymes, protein-S-isoprenylcysteine O-methyltransferase was encoded by two copies of the genes, whereas the remaining enzymes were each encoded by a single gene (Supplementary Table S3). Many genes involved in phenylalanine metabolism, benzoate degradation, toluene degradation, pyruvate metabolism, galactose metabolism, fatty acid degradation, and aminobenzoate degradation were identified in D. sinensis. For example, 24 key enzymes involved in the pyruvate metabolism pathway were identified (Supplementary Figure S4; Supplementary Table S4), with 37 genes encoding these enzymes, and 11 key enzymes involved in the phenylalanine metabolism pathway were identified, with 24 genes encoding these enzymes (Supplementary Figure S5; Supplementary Table S5).

The total number of genes encoding CAZymes among the two species of Daedaleopsis was more or less similar, as was the number of genes encoding each of the six classes of CAZymes (Table 2). Only a few differences could be concluded at the level of gene families encoding CAZymes (Figure 3). According to the number of genes belonging to different families of CAZymes, D. sinensis has a similar strategy to D. nitida for the utilization of woody substrates. From D. sinensis, 481 genes were assigned to the six classes of CAZymes, including 7 genes encoding carbohydrate-binding modules (CBMs), 33 carbohydrate esterases (CEs), 237 glycoside hydrolases (GHs), 73 glycosyltransferases (GTs), 18 polysaccharide lyases (PLs), and 113 auxiliary activities (AAs) (Table 2). Among the six classes, GHs were the dominant group (encoded by the most genes) and were mainly involved in the degradation of celluloses (GH1, GH3, GH5, GH6, GH7, GH9, GH12, GH51, and GH74), hemicelluloses (GH3, GH10, GH30, GH43, and GH51), pectins (GH28), chitins (GH18), and starches (GH31). The families each encoded by ten or more genes included AA2 (15), AA3 (29), AA5 (13), AA7 (20), AA9 (19), CE16 (13), GH16 (30), GH18 (21), GH28 (11), GH43 (10), GH5 (22), GH79 (12), and GT2 (20). The types and quantities of enzymes found in D. nitida were similar to those found in D. sinensis (Table 2; Supplementary Table S6).

CYPs play a variety of roles in fungi, including detoxification, heterobiomass degradation, and biosynthesis of secondary metabolites. Genes encoding CYPs comprise 4.3% of the CDS of the genome. After comparison with the CYP database, 524 genes encoding 39 CYPs were identified in D. sinensis (Table 2; Figure 4). In D. nitida, 40 CYPs encoded by 992 genes were identified. The most abundant CYP families in D. sinensis were CYP51, CYP620, CYP53, and CYP4.

A total of 44 gene clusters were predicted from D. sinensis. Of them, 23 encode terpene synthases (TS), 3 encode iterative type I polyketide synthases (T1PKS), 2 encode nonribosomal peptide synthetases (NRPS), 11 encode NRPS-like fragment, and 5 encode fungal unspecified ribosomally synthesized and post-translationally modified peptide product (fungal-RiPP-like) (Table 2). In its congeneric species D. nitida, 44 secondary metabolite clusters were predicted, of which 21 were TS, 4 were T1PKS, and the rest were the same as D. sinensis.

A total of 2,260 orthologous groups, including 1,402 single-copy genes, were identified from all 17 studied fungal species. The phylogenomic tree inferred from an alignment of the 1,402 single-copy orthologous genes with 530,910 characters of amino acid residues delimited the phylogenetic relationships among the 17 species with full bootstrap support (Figure 5).

Besides two species in Daedaleopsis, other white rot fungal species, viz. Trametes versicolor, Fomes fomentarius, Ganoderma sinense, Irpex rosettiformis, Pleurotus ostreatus, and P. chrysosporium, were also phylogenetically separated. Of the species in Daedaleopsis, D. sinensis and D. nitida, which occurred in a mean crown age of 28.77 Mya with a 95% highest posterior density of 21.22–36.67 Mya, had a close phylogenetic relationship. This is also the first time that the pairwise divergence time between these two species has been given. In the evolutionary process of the 17 sampled fungal species, gene family contraction occurred more commonly than gene family expansion (Figure 5). Regarding Daedaleopsis, 546 and 845 gene families had expanded in D. sinensis and D. nitida, respectively, whereas 1,348 and 573gene families had contracted in D. sinensis and D. nitida, respectively. Of these gene families, 46 (10 expanded and 36 contracted), and 59 (52 expanded and 7 contracted) in D. sinensis and D. nitida, respectively, experienced a rapid evolution.

Arrangement of the homologous genes or sequences between contigs >2.5 Mb in D. sinensis and contigs >0.6 Mb in D. nitida (Figure 6). Figure 6 shows only the connecting lines with similarity >5,000. Comparisons of homology allow for the study of evolutionary relationships between the two species.

Using the KOG database, 3,367 (27.71%) genes in D. sinensis were assigned to KOG categories (Figure 7). The most gene-rich KOG classifications were “general function prediction only,” “posttranslational modification, protein turnover, chaperones,” “secondary metabolites biosynthesis, transport, and catabolism,” and “signal transduction mechanisms.”

According to the GO database, 5,612 genes that accounted for 46.18% of the entire genome were distributed in the three functional categories of biological process, cellular components, and molecular function (Figure 8). Among biological processes, 17,835 genes were involved in “metabolic activities,” followed by “cellular process” (7,116 genes), “single-organism process” (2,627 genes), and “biological regulation” (1,617 genes) functions. Of the genes with products predicted to be involved in molecular function, 9,380, 2,407, and 329 were involved in binding, catalytic activity, and transporter activity, respectively. For cellular components, 3,224, 2,389, 1,996, and 1,029 genes were involved in the formation process of cell part, organelle, cell, and membrane, respectively.

In the two species of Daedaleopsis, the number of gene clusters involved in the synthesis of NRPS, NRPS-like, fungal-RiPP-like, and β-lactone was similar, and D. sinensis had the highest and the lowest numbers of gene clusters, which were involved in the syntheses of terpene and β-lactone, respectively (Table 2). In D. sinensis, there were 24 gene clusters encoding terpenes, which are key biosynthetic enzymes; 3 gene clusters encoding T1PKS, which are associated with the biosynthesis of polyketides; and 2 gene clusters encoding NRPSs, which are modular enzymes.

The genome sequencing and annotation of D. sinensis are crucial for its function and comparative genomics research. Following accurate species identification, this study presented the first whole genome sequencing of D. sinensis. Information on de novo genome assembly and characterization of D. sinensis is summarized in Table 1.

The most abundant CYP families in D. sinensis are CYP51, CYP620, CYP53, CYP4, CYP504, CYP505, CYP102, and CYP125. CYP51 is a structurally and functionally conserved fungal P450 family (Yu et al., 2020), and is widely distributed in different biological kingdoms, being found in animals, plants, fungi, yeast, lower eukaryotes, and bacteria (Yoshida et al., 2000; Lepesheva and Waterman, 2004). CYP51 is a key enzyme in biosterol synthesis (Yoshida et al., 2000) and can catalyze the 14α-methyl hydroxylation of the precursor of sterol (Waterman and Lepesheva, 2005; Lepesheva and Waterman, 2007). CYP51 has become an important target for cholesterol-lowering drugs, antifungal drugs, and herbicides (Lepesheva and Waterman, 2004). Genes in the CYP53 family are involved in degradation or detoxification of benzoate and its derivatives (Crešnar and Petrič, 2011; Yu et al., 2020). CYP505 genes are involved in ω-1 to ω-3 carbon hydroxylation of fatty acids (Nakayama et al., 1996; Bao et al., 2013). Search on biosynthesis of secondary metabolites showed that the CYP620 family is related to “carbohydrate metabolism-amino sugar and nucleotide sugar metabolism” and may play a role in terpenoid synthases (Yu et al., 2020). CYP504 family members are linked to the degradation of phenylacetate and its derivatives (Ferrer-Sevillano and Fernández-Cañón, 2007; Crešnar and Petrič, 2011).

According to the CAZyme database, GH genes involved in the degradation of glucoside bonds between sugars and sugars or between sugars and non-sugar groups were the most abundant CAZyme groups in D. sinensis, accounting for 49.27% of the total CAZyme gene sequences (Table 2). The second largest group with a proportion of 23.49% was the coenzyme family enzymes (AAs), which are involved in other CAZymes (some families are also involved in the degradation of lignin). Genes encoding GTs, which catalyze the biosynthesis of monosaccharides, disaccharides, oligosaccharides, and polysaccharides, accounted for 15.18% of the total CAZyme gene sequences in D. sinensis. CEs, which are responsible for sugar modification, and CBMs, which promote anchoring of lignocellulose-degrading enzymes to the corresponding substrate, accounted for 6.86% and 1.46%, respectively. The proportion of PLs, which are involved in the cleavage of polysaccharide aldehyde polyglycan chains, was 3.74%. In summary, the cumulative proportion of GH and PL genes responsible for the degradation of carbohydrates in D. sinensis was as high as 53.01%, and this reflects, to some extent, the ability of D. sinensis to degrade lignocellulose.

According to the KEGG database analysis of D. sinensis data, 3,831 genes are involved in nearly 120 metabolic pathways. These pathways mainly include carbon metabolism, fatty acid biosynthesis and degradation, amino and nucleotide sugar metabolism, terpenoid metabolism, starch and sucrose metabolism, metabolism of cofactors and vitamins, methane metabolism, polycyclic aromatic hydrocarbon degradation, pyruvate metabolism, phenylalanine metabolism, and carbohydrate metabolism. Among them, the metabolism of sugar-containing compounds includes the tricarboxylic acid cycle and the classic pentose phosphate pathway. In addition, fructose, mannose, and galactose metabolism pathways are predicted. Studies have shown that the degradation of lignin by white rot fungi generally belongs to secondary metabolism, whereas the degradation of cellulose and hemicellulose exists in both primary and secondary metabolic stages (Zhao et al., 2019). The KEGG results indicate that substrate and energy metabolism and trophic growth are among the major biological processes carried out by D. sinensis, which is consistent with the reality that catabolism and trophic growth are the predominant life activities of D. sinensis itself as a degrader. The genome data suggested that D. sinensis possesses several carbon metabolic pathways and glucose metabolic pathways related to the degradation of cellulose and hemicellulose. There are 35 genes in the genome involved in glycolysis or gluconeogenesis pathways, and 43 genes involved in starch and sucrose metabolism pathways. The starch and sucrose metabolic pathways are closely related to the degradation of cellulose. In the pathways, a total of 11 genes encoding cellulase were identified. Among them, eight genes encode β-glucosidase, two genes encode exoglucanase, and one gene encodes endoglucanase. These annotated cellulase genes are involved in the transformation of cellulose into _D-glucose (Supplementary Figure S1). Cellulose can be converted to cellodextrin by the action of endoglucanase, and then to cellobiose with the participation of β-glucosidase (or a combination of endoglucanase and exglucanase), which is hydrolyzed to _D-glucose by the action of β-glucosidase. The benzoic acid metabolism pathway is the main metabolic pathway in the process of lignin degradation, but this has been less studied in fungi and mainly reported in bacteria (Egland et al., 1997). In D. sinensis, the pathways associated with lignin degradation include those associated with the degradation of aromatic compounds, such as phenylalanine metabolism, benzoate degradation, toluene degradation, pyruvate metabolism, galactose metabolism, fatty acid degradation, and aminobenzoate degradation.