The dikaryotic G. lucidum strain GL0102 (the same strain with “Zhi 102” of Mycological Research Center, Fujian Agriculture and Forestry University) was maintained on Difco^™ Potato Dextrose Agar (PDA) at 4°C and stored at the Institute of Biological and Medical Engineering, Guangdong Academy of Sciences. Minimal medium plates containing glucose (20 g/l, MM), (NH4)₂SO₄ (1.5 g/l), K₂HPO₄ (1.0 g/l), MgSO₄ (0.3 g/l), vitamin B1 (0.5 mg/l), and agar (18 g/l) were prepared. Then, the glucose was replaced with an equal amount of 17 kinds of other carbon sources (20 g/l, Table 1) to make media with different carbon sources. Specifically, the AFRs (such as rice-straw and wheat-straw) were dried and crushed using a pulverizer, and were added to 500 ml triangular flasks together with other components of the media. One rotor was put in each triangular flask before sterilization. Then the triangular flasks were stirred on a magnetic stirrer for 20 s. The media (<50°C) were poured rapidly into the plates (9 cm in diameter), and stirred for 3 s after each pouring to ensure the homogeneity of the medium. Besides, to speed up the solidification and avoid settling of insoluble components, the plate cover was open and the media was blown by sterile air. The medium preparation ensured maximal contact of the mycelia with the carbon sources. GL0102 was then inoculated into the above culture media and PDA plates and cultured at 28°C for 5 days. For media with MM, lignin (LM) and cellulose (CM), 4 mycelial blocks were inoculated in each plate, and 60 plates in total; for media with other carbon sources, one mycelial block was inoculated in each plate, and 30 plates in total. For each sample, the mycelia (approximately 10–20 plates) were quickly scraped and mixed to produce a biological replicate, which was then frozen in liquid nitrogen and stored at −80°C. Three replicates were prepared for each of the treatments. Besides, GL0102 mycelia cultured on the following six carbon sources for 3 days and 7 days were also collected in the same way: glucose, lignin, cottonseed hull, rice-straw, bagasse, and PDA. The growth rates of GL0102 on different carbon sources were determined at 28°C based on the hyphal radical length.

RNA-Seq was performed from mycelia cultured on the above 19 different conditions for 5 days, each in triplicate. Extraction and quality control of the total RNA from each sample were conducted using methods previously reported (Wang et al., 2018), and RNA with RNA Integrity Number ≥ 7.5 was kept (Supplementary Table 1). cDNA library construction and sequencing were carried out following protocols of MGI sequencing platform and paired-end 150-bp reads were generated. Raw read quality was assessed using FastQC¹, and low-quality bases or reads were filtered out using Skewer (Jiang et al., 2014) with parameters: “-q 20 -Q 30 -l 50” (Supplementary Table 2). Expression levels (transcript per million, TPM) of genes in mycelia cultured on different carbon sources were calculated using HISAT2 and StringTie (Pertea et al., 2016). Two reference haploid genomes, GL0102_8 and GL0102_53 (two mating-compatible monokaryons from GL0102)² were downloaded from GPGD (Liao et al., 2021). Hierarchical clustering analysis of TPM values was performed using the “pheatmap” package³ in R. Differentially expressed genes (DEGs) were identified using the “DESeq2” package (Wang et al., 2010), using genome of GL0102_8 as the reference. DEGs were classified as genes that had log₂FoldChange ≥2 and value of p ≤1e-6.

Gene ontology (GO) enrichment analysis was carried out using clusterProfiler (Yu et al., 2012), and enrichment results with a value of p <1 × 10⁻³ were retained.

Ground powder (0.2 g) from mycelia in the presence of liquid nitrogen was extracted using 800 μl of 0.9% normal saline at 4°C for 10 min. After centrifuging (10,000 × g, 15 min, 4°C), the supernatant was used as the crude enzyme. The enzymatic activity of laccases was assayed by the ABTS (2, 2′-azino-bis-3-ethylbenzthiazoline-6-sulphonate) method (Johannes and Majcherczyk, 2000) using the following steps: 250 μl enzyme solution was added to the mixture of 500 μl acetate buffer (pH 4.5) and 250 μl ABTS (0.5 mm), then the absorbance at 420 nm was recorded within 3 min. One unit of enzyme activity was defined as the amount of laccase that oxidized 1 μm ABTS per minute.

Firstly, all the annotated proteins in GL0102_8 and GL0102_53 genomes were searched against the PFAM database (Pfam 32.0) using PfamScan (E value ≤1e-5).⁴ Genes with hits to PFAM ID PF00394.24, PF07731.16, and PF07732.17 were considered as the candidate laccase genes. Secondly, the genes were viewed and corrected using the Apollo browser (Dunn et al., 2019) as previously described (Wang et al., 2020). Thirdly, to rule out the false-positive results, all the laccase genes were reconfirmed using PfamScan. Finally, the whole set of laccase genes of GL0102_8 was searched against that of GL0102_53 using BlastP (Camacho et al., 2009). The allelic (one-to-one) relationship among these laccase genes was confirmed according to the amino acid (aa) similarity (identity ≥95%). For the other three individuals of G. lucidum, namely CGMCC5.0026 (Chen et al., 2012), P9 (single basidiospore; Liu et al., 2012), and Lingjian-2 (Tian et al., 2021), the laccase genes was identified following the same process above.

The number of aa, isoelectric point (pI), and molecular weight (MW) were computed using ProtParam.⁵ The presence of signal peptide cleavage sites was obtained using SignalP 5.0.⁶ The domains of laccases were identified using PfamScan.⁷ The duplication of laccase genes was identified by the two criteria: (1) the protein length of the shorter sequence covered ≥75% of the longer sequence and (2) the similarity of the two aligned sequences was ≥90%.

A multiple sequence alignment of laccase genes in GL0102_8 and GL0102_53 was performed on full-length protein sequences using clustal (Sievers et al., 2011). A maximum likelihood (ML) phylogenetic tree was constructed using RAxML (Stamatakis, 2014) with 1,000 bootstrap replicates. In addition, a ML phylogenetic tree containing laccase protein sequences in three additional G. lucidum individuals (CGMCC5.0026, P9, and Lingjian-2) was built following the same process.

Genetic variations were detected using MEGA7 (Kumar et al., 2016) after alignment using clustal, and the exhibition of SNPs were plotted using trackViewer (Ou and Zhu, 2019). To validate the SNP loci, the cDNA derived from mycelia cultured on RSM was used as PCR amplification template. The synthesizing of primers (Supplementary Table 3), cloning and sequencing of amplification products were conducted by Tsingke Biotechnology Co., Ltd.

RNA samples for sequencing were further synthesized into cDNA. Glyceraldehyde 3-phosphate dehydrogenase (gapdh) was used as a reference. Primers for the laccase genes (Supplementary Table 3) were designed and synthesized by Tsingke Biotechnology Co., Ltd. The qPCR reaction was performed using the Applied Biosystems ABI 7500 PCR System (ABI, United States) according to previously reported methods (Wang et al., 2017). The PCR amplification mixture contained 1 μl of cDNA, 10 μl of ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd), 0.4 μl of 10 μm forward and reverse primers, and 8.2 μl ddH₂O. The PCR reaction was performed with the initial denaturation step for 3 min at 95°C; 40 cycles of 3 s at 95°C and annealing at 60°C for 32 s, and a holding step for 30 s at 72°C. The melting curve (60–95°C) and agarose gel electrophoresis of qPCR products were used to check the specificity of each qPCR reaction. The standard curves were generated using a twofold dilution gradient of the cDNA. Amplification efficiencies (E = 10^-1/slope–1) and correlation coefficient (R²) values were calculated by standard curves. The relative gene expression was calculated with the 2^–ΔΔC_T method (Livak and Schmittgen, 2001). Parametric one-way analysis of variance (ANOVA) followed by Duncan’s test was used to calculate significant differences among different groups (p < 0.05).

The allelic expression bias was identified based on RNA-Seq read depths of SNP loci between alleles. The validation of allelic expression bias was conducted as follows: (1) primers were designed at flanking regions of SNP locus (Supplementary Table 3); (2) PCR amplification was conducted using cDNA derived from mycelia cultured on RSM; and (3) PCR products was cloned and then randomly selected and sequenced at Tsingke Biotechnology Co., Ltd.

Co-expression analysis was conducted by Multiscale Embedded Gene Co-expression Network Analysis (MEGENA) pipeline (Song and Zhang, 2015) as follows: (1) pairwise Pearson correlation coefficients of all the DEGs were calculated based on expression TPM values; (2) planar filtered network was set up by taking significant correlation pairs; and (3) multiple correspondence analysis was employed to cluster genes into modules.

The 3 kb upstream of the start codon of laccase was extracted as the promoter region. Promoter sequences were aligned and analyzed by MEGA7. The binding sites matrix of fungi were downloaded from JASPAR⁸, and CiiiDER⁹ were used to predict the possible binding sites.

The G. lucidum strain GL0102 was cultured on a total of 19 kinds of media with different carbon sources as solo carbon source (Table 1), respectively, and the growth rate was calculated accordingly. Mycelia cultured on media with glucose (MM), lignin (LM) or cellulose (CM) as carbon sources were sparse and had much lower growth rate (0.17 ± 0.004 to 0.37 ± 0.008 cm/d; Figure 1A; Supplementary Figure 1), while mycelia cultured on media with AFRs as carbon sources were lush and had higher growth rate, especially on wheat-straw (WSM), corncob (COM), pteridophyte (PTM), and chestnut leaf (CLM; 0.76 ± 0.015 to 0.80 ± 0.012 cm/d; Figure 1A; Supplementary Figure 1). This indicates that single-component carbon sources are possibly not comparable to complex carbon sources (AFRs) for mycelia production of G. lucidum.

To uncover the key genes involved in mycelia growth on different carbon sources, we compared the gene expression data obtained from the mycelia cultured on 18 carbon sources with that on glucose (MM), respectively. Among all the pairwise comparisons, 7 (XM vs. MM) to 1,865 (SNM vs. MM) genes were downregulated, while 16 (XM vs. MM) to 1,006 (RSM vs. MM) were upregulated (Figure 1B). According to their expression profile, all the DEGs could be divided into three main clusters (Supplementary Figure 2). Genes in cluster 1 showed relatively high expression on MM, XM, CM, and PDA; genes in cluster 2 showed relatively high expression on LM; and genes in cluster 3 showed relatively high expression on carbon sources derived from AFRs. Besides, the genes of these three clusters all showed enrichment in carbohydrate degradation-related functions, such as oxidoreductase activity, hydrolase activity, and carbohydrate metabolic process (Figure 1C).

Notably, the laccase gene family, which is closely related with lignin degradation, showed significant expression dynamics when the mycelia was cultured on different carbon sources. Thus, the laccase activity was tested. The mycelia cultured on single-component carbon source showed quite low laccase activity, such as 26 ± 7.3 U/mg on MM, 29 ± 13.1 U/mg on XM, and 32 ± 0.78 U/mg on CM (Figure 1A), while laccase activity of mycelia cultured on AFRs (except PTM) was significantly induced, showing values greater than 281 ± 35.4 U/mg (Figure 1A; Supplementary Figure 1). Specially, mycelia cultured on CSM showed laccase activity of 2,410 ± 122.4 U/mg which was much higher than those on other carbon sources (Figure 1A). Given that there are multiple laccase members in the G. lucidum genome, and the characteristics of each laccase member and their contribution to laccase activity may be different (Wang et al., 2015b; Ha et al., 2017; Kumar et al., 2017), comprehensive analysis regarding G. lucidum laccase was then performed.

Totally, 15 and 16 laccase genes were identified in GL0102_53 and GL0102_8, respectively. According to the sequence similarity, two laccase genes in GL0102_8 were speculated to derived from duplication event. Hence, laccase genes in the two genomes were divided into 15 categories (Gllac1-Gllac15), and the two duplicated laccase genes in GL0102_8 were separately named as H8lac7_1 and H8lac7_2 (Supplementary Table 4). The laccase genes located on seven chromosomes of GL0102_8 and GL0102_53 and emerged the most frequently on chromosome 6 (Supplementary Table 4). Paired allelic laccase genes showed highly similar sequence features in protein length, pI, and MW. For example, H8lac14 and H53lac14 had minimal protein length of 460 aa, and H8lac7_1 and H53lac7 had maximal protein length of 667 aa (Supplementary Table 4). All the laccase genes were predicted to encode secreted proteins.

A ML phylogenetic tree of laccase genes of GL0102 was constructed, and all paired allelic laccase genes clustered with each other (Figure 2A). It was reported that G. lucidum strain CGMCC5.0026 and P9 had 13 (Chen et al., 2012) and 16 (Liu et al., 2012) laccase genes, while 15 and 17 laccase genes were identified, respectively, after the manual correction in this study (Gllac5 and Gllac6 each had one more copy in P9; Supplementary Table 5). Besides, 15 laccase genes were identified in G. lucidum strain Lingjian-2 (Supplementary Table 5). In the phylogenetic tree constructed with all the above laccase genes, each of the 15 laccase categories from different strains/haploids were clustered together (Supplementary Figure 3), indicating that laccase genes were relatively conserved in G. lucidum.

As shown in Figure 2B, paired allelic laccases presented highly conserved gene structures and possessed relatively large number of introns, for example, each of Gllac8, Gllac9, and Gllac10 had 12 introns. Each laccase contained three conserved domains: Cu-oxidase (PF00394.24), Cu-oxidase_2 (PF07731.16), and Cu-oxidase_3 (PF07732.17; Figure 2C).

Genetic variations were detected between paired allelic laccase genes, and SNPs were identified as the most abundant variations (Supplementary Table 6). Among all the allelic laccase pairs, 1 to 48 SNPs causing 0 to 16 amino acid changes were observed in coding regions, while none of the amino acid changes occurred in laccase signature sequence regions (Supplementary Table 6; Supplementary Figure 4). In general, the SNP rate was higher in untranslated regions (UTR) and introns. For example, a total of 113 SNPs and 10 insertions/deletions were identified in genic regions between Gllac2 alleles, of which 42 SNPs existed in coding regions (Figure 2D). Randomly selected SNP loci between H53lac2 and H8lac2 were confirmed by Sanger sequence (Figure 2D).

The expression pattern of laccase genes on different carbon sources was further analyzed based on the RNA-Seq data using GL0102_8 genome as reference. Nine laccase genes (Gllac1, Gllac4, Gllac6, Gllac8, Gllac9, Gllac10, Gllac11, Gllac15, and H8lac7_2) showed quite low expression level (barely detectable, TPM value <10 under most conditions) under all the testing conditions (except that H8lac7_2 showed relatively high expression level on LM), while the other seven laccase genes (Gllac2, Gllac3, Gllac5, Gllac7, Gllac12, Gllac13 and Gllac14) showed high expression level especially on complex carbon sources (TPM value >100 under certain conditions; Figure 3A). Notably, Gllac2 and Gllac7 were the top two expressed laccase genes on complex carbon sources. Differential expression analysis was conducted on all carbon sources pairwisely. In all comparisons, 1 to 10 laccases were found differentially expressed. And the maximal number of differentially expressed laccase genes were found in the comparisons of XM vs. SNM and PDA vs. SNM (Supplementary Figure 5).

To test if the laccase expression was related with different growth stages, GL0102 mycelia cultured for 3, 5, and 7 days on MM, LM, CHM, RSM, BM, and PDA were collected, respectively. And the transcription levels of Gllac2, Gllac4, Gllac7, Gllac9, Gllac10, and Gllac12 were analyzed by qPCR after we evaluated multiple parameters of the qPCR reaction and confirmed their acceptability (Supplementary Figure 6). Overall, the expression levels of Gllac2 and Gllac7 were highest in mycelia cultured for 3 days (on LM, RSM, BM and CHM) and decreased gradually with the increase of culture time (Supplementary Figure 7). Specifically, compared to 5d-mycelia cultured on MM, Gllac2 showed 5-fold (on CLM) to 1,149-fold (on PNM) upregulation and Gllac7 showed 4-fold (on PTM) to 94-fold (on RSM) upregulation on complex carbon sources (Figures 3B,C). Gllac4, Gllac9 and Gllac10 maintained low expression levels from 3 to 7 days (Supplementary Figure 7).

Further analysis showed that the two nuclei in G. lucidum contributed unevenly to the total expression of laccase, in other words, allelic expression bias was observed in laccase. On all the 19 carbon sources except XM, the expression of H53lac7 was consistently higher than that of H8lac7_1 (value of p < 0.05), and the read count ratio of H53lac7 ranged from 51.71 to 76.54% with an average of 62.6% (Figure 3D; Supplementary Figure 8). To confirm this allelic bias case, partial sequences of Gllac7 (311 bp), which contains SNP loci, were amplified using cDNA derived from RSM-cultured mycelia (Figure 3E). After cloning and sequencing, 19 out of 33 clones were consistent with H53lac7, while 14 clones were consistent with H8lac7_1 (Figure 3F), confirming the existing of allelic expression bias. To analyze the possible causes of allelic expression bias, we analyzed the promoter regions of H8lac7_1 and H53lac7. Totally, 47 SNPs and three insertions/deletions with length of 5, 8 and 9 bp were identified. These variants may cause differences in binding sites of transcription factors. For example, zinc fingers, a type of important regulators, have many binding sites in the promoter of Gllac7, while differences of binding sites were observed between promoters of H53lac7 and H8lac7_1 (Figure 4A).

A co-expression network was built based on expression levels of DEGs (in GL0102_8) using MEGENA pipeline (Song and Zhang, 2015). In the co-expression network, Gllac7 (H8lac7_1), H8lac3, H8lac12, and H8lac13 were clustered into a shared module (Figure 4B). This module contained 299 genes, among which five were TFs. Overall, there were 68 TFs in the DEGs and most TFs were downregulated on complex carbon sources. Moreover, Pearson correlation analysis between Gllac7 and these TFs indicated that 35 TFs were significantly correlated with Gllac7, of which six TFs showed positive correlation, while the other 29 TFs showed negative correlation with Gllac7 (Figure 4C). Most of these TFs belong to zinc finger family such as C2H2 (7/35) and Zn2C6 (12/35).