De novo transcriptomic sequencing and reads assembly. In order to compare the gene expression differences between the NZF and XZF cultured fruit bodies, we constructed a cDNA library, and we sequenced and assembled the transcriptome. As shown in Table 2, the total number of reads generated by s NZF and XZF were 49,596,004 and 54,986,028, respectively. By removing the reads with lengths less than 50 bp and average sequence quality less than Q20, the final clean read numbers were 46,582,176 and 51,766,838, respectively. The Q20 and Q30 values of each sample were higher than 97% and 93%, respectively, indicative of high-quality sequencing data that could be analyzed in the next step.
Table 2
Statistical results Sample: sample name; Reads No: total number of Reads; Clean Reads No: high-quality serial reads; N (%): percentage of fuzzy bases; Q20(%): the percentage of bases whose base recognition accuracy was more than 99%; Q30(%): the percentage of bases whose base recognition accuracy was above 99.9%; Clean Reads %: the percentage of high-quality sequence reads compared to the total sequenced reads
Sample
|
Reads No
|
Clean Reads No
|
N(%)
|
Q20(%)
|
Q30(%)
|
Clean Reads(%)
|
NZF1
|
49,596,004
|
46,582,176
|
0.000544
|
97.62
|
93.92
|
93.92
|
XZF
|
54,986,028
|
51,766,838
|
0.00058
|
97.77
|
94.27
|
94.14
|
Function annotation classification. All single genes were compared to sequences from the NR, GO, KEGG, Pfam, eggNOG, and Swissprot databases using the BLASTx algorithm (E-value ≤ 1.0×10− 5), and annotations were obtained (Table 3). A total of 18,391 genes were effectively annotated, of which 11,435 were successfully annotated using the NR database, accounting for 62.18% of the unigene number, followed by Swissprot (11,155, 60.65%), eggNOG (10,539, 57.31%), Pfam (8,136, 44.24%), GO (6,947, 37.77%), and KEGG (5,994, 32.59%). The species distribution of the annotations was analyzed (Fig. 2). In the predicted protein homology distribution, most of the proteins (37.88%) had a 60–80% homology with proteins from other fungi in the NR database, and only a few (9.95%) had < 40% homology. According to the NR database, the gene sequences of A. cinnamomea had the largest number of matches with Fibropolia radiosa (14.54%) and Laetiporus subphureus 93 − 53 (13.17%), followed by Daedalea quercina L − 15889 (5.53%), Gelatoporia subpopura B (4.78%).
Table 3
Summary of annotation results. Database: type of database; Number: the number of unigenes successfully annotated in the respective database; Percentage (%): the proportion of unigenes successfully annotated in the database to the total unigenes; In all database: The number of unigenes annotated in all databases
Database
|
Number
|
Percentage (%)
|
NR
|
11,435
|
62.18
|
GO
|
6,947
|
37.77
|
KEGG
|
5,994
|
32.59
|
Pfam
|
8,136
|
44.24
|
eggNOG
|
10,539
|
57.31
|
Swissprot
|
11,155
|
60.65
|
In all database
|
3,661
|
19.91
|
Differentially expressed genes (DEGs). The differentially expressed genes were assessed using DESeq and volcano mapping with the R language ggplots2 package (Fig. 3 and Table S1). 592 DEGs were identified between NZF and XZF (170 were up-regulated, 422 were down-regulated). Among the up-regulated genes were two genes encoding cytochrome P450s (AC263_c0_g1, AC3065_c0_g1), one terpenoid synthase (AC3097_c0_g1), and one GST (AC3305_c0_g1). Among the down-regulated genes, 11 genes were found encoding for cytochrome P450s (AC49_c4_g1, AC218_c0_g1, AC245_c2ug1, AC363_c0_g1, AC588_c0_g1, AC1102_c1ug3, AC1414_c0_g1, AC1641_c0_g3, AC1685_c0_g1, AC2006_c0_g1, AC2915_c0_g1), 1 TPS (AC2874_c0_g1) and 1 acetyl CoA synthase (AC57_c0_g1). These findings indicate that several cytochrome P450 genes play an important role in triterpene synthesis.
GO functional enrichment analysis of DEGs. To explore the function of NZF and XZF unigenes, we performed GO enrichment analysis (Fig. 4 and Table S2). According to the international standardized gene function classification system25, a 8,679 unigenes were divided into 3 main functions (biological process, cell component, and molecular function) and 30 functional subgroups. 2,010 DEGs were enriched in the molecular function category (487 were up-regulated and 1,523 were down-regulated), and the dominant subgroups were "structural construction of ribosome" and "structural molecule activity". 4,573 DEGs were enriched in the biological process category, (1,037 were up-regulated, and 3,536 were down-regulated). The dominant subgroups were "amide biological process" and "translation." 2,096 DEGs were enriched in the cellular component category (301 were up-regulated, and 1,795 were down-regulated). Notably, except for one up-regulated gene enriched in the amide biological process, 261 enriched in the six dominant subgroups were down-regulated genes. And their expression in NZF was much higher than that in XZF. Most were genes encoding for ribosomal proteins, indicating that they play an important role in protein biosynthesis in the A. cinnamomea fruit body.
KEGG enrichment analysis of the DEGs. To better understand the biological function of the unigenes, we compared 4,959 unigene sequences with the reference KEGG database pathways. A total of 216 DEGs were assigned to 71 KEGG pathways (Fig. 5 and Table S3). In addition, the unigenes were classified into six pathway categories: Genetic Information Processing, Metabolism,、Environmental Information Processing, Cellular Processes, Organismal Systems, and Human Disease. Many unigenes in the A. cinnamomea fruiting bodies were enriched in Metabolism, including Carbohydrate metabolism, Lipid metabolism, Amino acid metabolism, Glycan biosynthesis and metabolism, Xenobiotics biodegradation and metabolism, and Metabolism of cofactors and vitamins. In the top 10 pathways with significant enrichment, 28 differentially expressed genes (3 up-regulated and 25 down-regulated) had a much higher expression in XZF, except for one gene encoding for endoglucanase. These results provide important insights for further research into the development and function pathways that are specific to NZF and XZF.
Differentially expressed terpenoid biosynthesis genes between NZF and XZF. The FPKM values were for gene expression quantification, and 44 genes were found to be differentially expressed between NZF and XZF (Fig. 6 and Table S4). Among the terpenoid backbone biosynthesis genes, ACS1, ACS2, and ACS5 were significantly highly expressed in NZF. Among them, ACS1 expression in NZF was two-fold higher compared to XZF. On the other hand, ACS3 and ACS4 expression was higher in XZF; however, the expression difference was not significant. The OSC, FPS, IPP, TPS, MPD, HMGCS, and XI-FIT genes were all highly expressed in XZF; however, the difference in expression was not significant. All 14 genes encoding glycosyltransferases (GT) were highly expressed in XZF; however, GT1, GT2, GT3, GT5, GT6, GT7, GT9, GT11, GT12, and GT13 exhibited non-significant differences in expression between NZF and XZF, The difference in expression of GT4, GT8, GT9, GT10, and GST1 was significant, and the expression of GT10 was 1.5-fold higher in XZF than that in NZF. CYP4501, CYP4502, CYP4503, CYP4504, CYP4507, CYP4512, and CYP4516 encoding cytochrome P450s were highly expressed in XZF, albeit with non-significant differences in expression. CYP61, CYP51, CYP4505, CYP4506, CYP4508 CYP4509, CYP4510, CYP4511, CYP4513, CYP4514, CYP4515 were significantly differentially expressed. CYP61 and CYP51 exhibited a two-fold higher expression in XZF compared to NZF, while CYP4506 and CYP4515 exhibited a five-fold higher expression in NZF compared to XZF. The above data indicate that the culture substrate influences the synthesis of terpene metabolites and suggests that C. camphora is viable as an alternative host of A. cinnamomea.
Validation of DEGs by qRT-PCR. To further validate the RNA-Seq data accuracy, we randomly selected 9 genes of the triterpene synthesis pathway and performed qRT-PCR to validate their expression in NZF and XZF. qRT-PCR results (Fig. 7) showed confirmed the differential expression patterns of these genes. With NZF as the control, the expression of IPS, CYP51, 2,3-OSC, TPS, and GT1, encoding the terpenoid backbone biosynthesis genes, were up-regulated in XZF, while genes encoding for CYP450, ACS, SQS, and GT2 were down-regulated, respectively. The above data further indicated that the culture substrate affected terpene metabolites synthesis, with their expression between NZF and XZF being highly variable.