微生物学报  2023, Vol. 63 Issue (10): 4016-4033   DOI: 10.13343/j.cnki.wsxb.20230150.
http://dx.doi.org/10.13343/j.cnki.wsxb.20230150
中国科学院微生物研究所,中国微生物学会

文章信息

肖冬来, 马璐, 杨驰, 刘晓瑜, 林辉, 江晓凌. 2023
XIAO Donglai, MA Lu, YANG Chi, LIU Xiaoyu, LIN Hui, JIANG Xiaoling.
基于转录组和加权基因共表达网络的广叶绣球菌木质纤维素降解特性分析
Characterization of lignocellulose degradation by Sparassis latifolia based on transcriptomics and weighted gene co-expression network analysis
微生物学报, 63(10): 4016-4033
Acta Microbiologica Sinica, 63(10): 4016-4033

文章历史

收稿日期:2023-03-09
网络出版日期:2023-07-17
基于转录组和加权基因共表达网络的广叶绣球菌木质纤维素降解特性分析
肖冬来 , 马璐 , 杨驰 , 刘晓瑜 , 林辉 , 江晓凌     
福建省农业科学院食用菌研究所, 福建 福州 350014
摘要[目的] 分析广叶绣球菌(Sparassis latifolia)在不同木质纤维素诱导条件下基因表达差异,为广叶绣球菌木质纤维素降解关键基因和分子机制研究提供参考。[方法] 以松木、杉木、甘蔗渣和天然堆积发酵后的杉木和发酵后的甘蔗渣为碳源,在液体培养条件下培养诱导广叶绣球菌,对其转录组进行测序研究,并对不同木质纤维素诱导样本进行加权基因共表达网络分析(weighted gene co-expression network analysis, WGCNA)。[结果] 杉木培养与松木培养比较组差异表达基因最少(20个),蔗渣培养与松木培养比较组差异表达基因最多(486个)。基因本体(gene ontology, GO)富集分析结果表明,差异表达基因主要涉及氧化还原酶活性、单加氧酶活性和铁离子结合活性等,京都基因和基因组百科全书(Kyoto encyclopedia of genes and genomes, KEGG)通路富集分析结果表明,差异表达基因主要涉及戊糖和葡萄糖醛酸转换、甲烷代谢和乙醛酸盐和二羧酸盐代谢等通路。发酵甘蔗渣为碳源培养时,纤维素和半纤维素降解相关的糖苷水解酶基因表达量总体上较高,而未发酵的松木、杉木和甘蔗渣为碳源培养时木质素降解或修饰相关的碳水化合物辅助酶基因表达量总体上较高。利用WGCNA共鉴定出10个共表达模块,其中green模块与未发酵蔗渣诱导显著正相关,blue模块与发酵甘蔗渣诱导显著正相关,magenta和turquoise模块与发酵杉木诱导显著正相关。GO富集分析结果表明,turquoise模块内基因显著富集到尿素跨膜转运子活性、甲基转移酶活性和单加酶活性等,blue模块基因显著富集到水解酶活性和β-甘露糖苷酶活性。KEGG通路富集分析结果表明,blue模块内基因显著富集的通路有半乳糖代谢、果糖和甘露糖代谢、苯丙氨酸代谢、精氨酸和脯氨酸代谢等。通过构建互作网络图挖掘到12个核心基因,其可能参与了基质降解及相关基因的表达调控。[结论] 不同木质纤维素类型显著影响了广叶绣球菌木质纤维素降解基因的差异表达轮廓,这种差异反映了广叶绣球菌对不同木质纤维素特异的降解策略。
关键词广叶绣球菌    转录组分析    木质纤维素降解    加权基因共表达网络分析    
Characterization of lignocellulose degradation by Sparassis latifolia based on transcriptomics and weighted gene co-expression network analysis
XIAO Donglai , MA Lu , YANG Chi , LIU Xiaoyu , LIN Hui , JIANG Xiaoling     
Institute of Edible Fungi, Fujian Academy of Agricultural Sciences, Fuzhou 350014, Fujian, China
Abstract: [Objective] To explore gene expression profiles of Sparassis latifolia cultivated with different lignocellulose materials and provide a reference for mining the key genes and deciphering the mechanism of lignocellulose degradation by S. latifolia. [Methods] We examined the transcriptomes of S. latifolia cultivated in the liquid medium with pine, spruce, bagasse, fermented spruce, or fermented bagasse as the sole carbon source. Weighted gene co-expression network analysis (WGCNA) was performed on the gene expression of samples cultivated with different lignocellulose materials. [Results] There were only 20 differentially expressed genes (DEGs) between the spruce and pine groups, and 486 DEGs (the highest number) between bagasse and pine groups. Gene ontology (GO) enrichment analysis showed that the DEGs were mainly involved in oxidoreductase activity, monooxygenase activity, and iron ion binding. Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis showed that the DEGs were mainly involved in pentose and glucuronate interconversions, methane metabolism, and glyoxylate and dicarboxylate metabolism. The expression of genes encoding glycoside hydrolases associated with cellulose or hemicellulose degradation was generally higher when the strain was cultured with fermented bagasse as the sole carbon source. The expression of genes encoding carbohydrate-active enzymes for lignin degradation or modification was generally higher when the strain was cultured with pine, spruce, or bagasse as the sole carbon source. A total of 10 co-expression modules were identified by WGCNA. Significant positive correlations existed between the green module and bagasse, between the blue module and fermented bagasse, and between magenta and turquoise modules and fermented spruce. GO enrichment analysis showed that the genes in the turquoise module were enriched in urea transmembrane transporter activity, methyltransferase activity, and monooxygenase activity, while those in the blue module in hydrolase activity and β-mannosidase activity. KEGG pathway enrichment analysis showed that the genes in the blue module were enriched in galactose metabolism, fructose and mannose metabolism, phenylalanine metabolism, arginine and proline metabolism, etc. From the interaction network, 12 hub genes were obtained, which may be involved in lignocellulose degradation or gene expression regulation. [Conclusion] Different lignocellulose materials significantly affected the gene expression profiles of S. latifolia, which may imply the specific degradation strategies against different lignocellulose materials.
Keywords: Sparassis latifolia    transcriptomic analysis    lignocellulose degradation    weighted gene co-expression network    

褐腐菌在木质纤维素降解过程中主要依靠木质纤维素降解酶系,同时进行一些非酶促反应如芬顿反应(Fenton’s reaction)和哈伯-韦斯反应(Haber-Weiss reaction)等,这些非酶反应有助于破坏植物细胞壁木质素结构暴露出纤维素、半纤维素,在协助基质降解过程中发挥重要作用[1-2]。木腐菌可利用的木质纤维素范围存在较大差异,某些物种只在特定的植物上出现,而有些则具有广泛的基质利用范围[3-4]。不同基质的木质纤维素组成、化合物种类及含量差异影响着这些真菌对基质的降解[5]。木质纤维素类型对木腐真菌降解酶系的诱导表达及RNA水平的编辑存在基质特异性[6]。相对于简单碳源(葡萄糖),木质纤维素类复杂碳源可诱导碳水化合物活性酶(carbohydrate active enzymes, CAZymes)基因及与芬顿反应相关的萜烯合成酶和细胞色素P450等基因的上调表达[7]。在蛋白质组水平上,CAZymes的表达同样受到不同木质纤维素种类的影响[8-9]

广叶绣球菌(Sparassis latifolia)主要生长在针叶树上,如落叶松(Larix kaempferi)、赤松(Pinus densiflora)和红松(Pinus koraiensis)等[10]。食用菌栽培基质很大部分来源于农林废弃物,如木屑、棉籽壳、玉米芯和甘蔗渣等。目前,国内绣球菌栽培主料为松木屑。课题组前期的栽培基质筛选结果表明绣球菌在单一杉木、甘蔗渣以及自然发酵后的杉木、自然发酵后的甘蔗渣上均可正常生长。

目前,有关绣球菌对木质纤维素降解的分子机制研究还较少[11-12]。本研究利用转录组测序分析了5种不同类型木质纤维素(松木屑、杉木屑、未发酵甘蔗渣和自然条件下建堆发酵6个月左右的杉木屑及甘蔗渣)在液体培养条件下诱导广叶绣球菌的基因表达差异。进一步结合加权基因共表达相关网络分析(weighted gene co-expression network analysis, WGCNA),以挖掘该菌与不同木质纤维素类型相关的关键模块和基因。研究结果将为广叶绣球菌特异性木质纤维素降解功能基因和降解机制的研究提供基础。

1 材料与方法 1.1 试验菌株和培养条件

广叶绣球菌S. latifolia (闽认菌2013005)为福建省农业科学院食用菌研究所保藏。菌株培养、活化及样品制备参考杨驰等[12]的方法,并稍作修改。利用打孔器(直径2 mm)将经马铃薯葡萄糖琼脂培养基(potato dextrose agar, PDA)活化的菌块接种于100 mL种子培养基(土豆20%,葡萄糖2%,鱼粉蛋白胨0.2%,pH值5.5),25 ℃、避光振荡培养9 d (150 r/min)。

1.2 木质纤维素样品及成分测定

未发酵新鲜松木屑(标记为Pine)、杉木屑(Spru)、甘蔗渣(Bag)以及室外建堆自然发酵6个月左右的杉木屑(Spru-FD)和建堆自然发酵6个月左右的甘蔗渣(Bag-FD)均来源于福建省南平市顺昌县食用菌生产企业,室外建堆发酵方法参考文献[13]。将以上木屑和甘蔗渣烘干后分别粉碎,过9目筛后收集截留部分(长度约2 mm)作为样品备用;称取0.50 g样品测定灰分含量,并依次测定中性洗涤后纤维(neutral detergent fiber, NDF)、酸性洗涤后纤维(acid detergent fiber, ADF)和酸性洗涤后木质素(acid detergent lignin, ADL)含量[14],按公式分别折算半纤维素、纤维素和木质素的含量。

利用SPSS软件,采用Waller-Duncana法进行多组样本间组分含量差异显著性分析。

1.3 转录组试验样品制备

转录组试验样品制备依据课题组前期关于差异基因表达的研究结果加以设计[11-12],将不同木质纤维素样品细粉经126 ℃灭菌120 min后,分别称取0.5 g加入含100 mL无菌水的250 mL三角瓶中作为诱导培养基。将1.1中振荡培养9 d的菌丝体分别加入到不同的诱导培养基中继续培养5 d (25 ℃, 150 r/min),过滤去除不同碳源细粉后获得菌丝体,液氮速冻后送上海美吉生物医药科技有限公司进行转录组测序分析。

1.4 转录组测序分析

使用MJzol总RNA提取试剂盒(上海美吉生物医药科技有限公司)提取不同诱导处理的绣球菌总RNA。RNA质检合格后构建文库上机测序。使用fastp[15]对测序原始数据进行过滤,获得高质量测序数据。与绣球菌参考基因组[16]比对后,利用RSEM[17]软件计算不同试验处理样品中转录本丰度。利用DESeq2[18]软件分析不同基因表达量差异倍数(fold change, FC),以|log2 FC|≥1且经过校正的P值(P-adjust) < 0.05为标准筛选获得发酵杉木/松木、杉木/松木、蔗渣/松木、发酵蔗渣/松木、发酵蔗渣/未发酵蔗渣和发酵杉木/未发酵杉木比较组间的差异表达基因。分别采用Goatools和KOBAS[19]软件对差异表达基因进行基因本体(gene ontology, GO)功能富集和京都基因和基因组百科全书(Kyoto encyclopedia of genes and genomes, KEGG)通路富集分析,P-adjust < 0.05时为显著富集。通过美吉生物单因子相关性分析云工具(https://cloud.majorbio.com/page/tools/)分析不同基质木质纤维素组分含量与差异基因表达相关性(Spearman相关系数≥0.85,且P<0.01)。

1.5 加权基因共表达网络分析

WGCNA分析利用美吉生物云平台(https://cloud.majorbio.com/)完成。模块内最少基因数设置为30,相似模块合并阈值设置为0.25。采用Pearson相关性检验分析模块与性状关系,相关系数≥0.8,且P<0.05为显著正相关。网络生成后,选取模块内连通性前30的关键节点使用Cytoscape 3.9.1及CytoNCA插件绘制网络图。

2 结果与分析 2.1 不同生物质木质纤维素组分分析

检测分析了5种不同生物质原料中的主要成分,结果见表 1。从表 1可以看出,甘蔗渣(Bag)和发酵甘蔗渣(Bag-FD)纤维素、半纤维素含量高于松木(Pine)、杉木(Spru)和发酵杉木(Spru-FD),而木质素含量要低于松木和杉木。甘蔗渣和发酵甘蔗渣间、松木和杉木间纤维素含量差异不显著,松木和发酵杉木间半纤维素含量差异不显著,杉木和发酵杉木间木质素含量差异不显著。

表 1. 五种不同生物质原料化学组成分析 Table 1. The chemical properties of five kinds of biomass materials
Sample Cellulose (%) Hemicellulose (%) Lignin (%) Ash (%)
Pine 48.75±0.91bc 10.51±0.21c 29.00±0.43b 5.17±0.15a
Spru 49.71±1.23b 7.37±0.63d 33.26±0.94a 2.36±0.04d
Spru-FD 46.19±0.42c 9.88±1.43c 33.90±1.04a 4.53±0.10b
Bag-FD 57.04±2.15a 21.44±0.07a 14.97±0.17c 3.39±0.11c
Bag 58.65±1.73a 18.19±0.39b 12.64±1.03d 1.38±0.04e
Different lowercase letters in the same column indicate significant differences (P<0.05).

2.2 不同生物质诱导差异表达基因的鉴定

通过转录组测序分析了5种不同类型基质(松木、杉木、未发酵甘蔗渣和自然条件下建堆发酵的杉木屑及甘蔗渣)分别作为单一碳源时,在液体培养条件下诱导广叶绣球菌的基因表达差异。以|log2 FC|≥1且P-adjust < 0.05为标准,获得5种基质碳源培养下,不同2组基质间差异表达基因(differentially expressed genes, DEGs),结果如图 1所示。从图 1A看出,发酵杉木培养与松木培养相比(Spru-FD/Pine)获得DEGs 219个,杉木培养与松木培养相比(Spru/Pine)获得DEGs 20个,蔗渣培养与松木培养相比(Bag/Pine)获得DEGs 486个,发酵蔗渣培养与松木培养相比(Bag-FD/Pine)获得DEGs 357个,发酵蔗渣培养与未发酵蔗渣培养相比(Bag-FD/Bag)获得DEGs 463个,发酵杉木培养与未发酵杉木培养相比(Spru-FD/Spru)获得DEGs 383个。杉木培养与松木培养相比DEGs最少,其中上调表达5个,下调表达15个。蔗渣培养与松木培养相比DEGs最多,其中上调表达399个,下调表达87个。从图 1B看出,以松木培养的转录组为对照,杉木、未发酵甘蔗渣、发酵杉木及发酵甘蔗渣培养的转录组间共有的差异基因数为3个(EVM0003544、EVM0007291和EVM0011928),分别编码乙酰辅酶A合成酶类似蛋白(acetyl-CoA synthetase-like protein)、锌型乙醇脱氢酶类似蛋白(zinc-type alcohol dehydrogenase-like protein)和芳基醇脱氢酶(aryl-alcohol dehydrogenase)。

图 1 差异表达基因统计图(A)和韦恩图(B) Figure 1 Statistical histogram (A) and Venn diagram (B) of differentially expressed genes.

通过5种不同生物质木质纤维素组成含量与差异表达基因相关性分析(Spearman相关系数≥0.85,且P<0.01),共获得与木质素含量正相关的基因2个,与纤维素含量正相关的基因8个,与半纤维素含量正相关的基因17个。如图 2所示,与木质素含量正相关的亚分支3 (subcluster 3)中的2个基因功能未知,编码假定的功能未知蛋白;与半纤维素含量正相关的亚分支2 (subcluster 2)包含3个糖转运蛋白(sugar transporter),1个MFS通用底物转运蛋白(MFSgeneral substrate transporter)以及2个短链脱氢酶(short chain dehydrogenase)等蛋白编码基因。与纤维素含量正相关的亚分支1 (subcluster 1)包含锌型乙醇脱氢酶类似蛋白、MFS通用底物转运蛋白以及糖苷水解酶家族115 (glycoside hydrolase family 115)等蛋白编码基因。其中假定蛋白EVM0007776和糖转运蛋白EVM0006328编码基因分别与木质素和半纤维素含量相关性较高,相关系数(correlation coefficient)大于0.9。

图 2 木质纤维素组分显著相关差异基因热图分析 Figure 2 Heat map analysis of DEGs significantly associated with differences in lignocellulosic components.

2.3 差异表达基因的GO和KEGG富集分析

利用Goatools和KOBAS软件对不同比较组间DEGs进行了功能富集分析。其中杉木培养与松木培养比较组未显著富集到相应的GO条目,说明可能由于差异基因数量较少,相关基因表达的改变导致调控途径原有功能的变化不显著,其余比较组的GO富集结果见图 3A。从图 3A可以看出,发酵杉木培养与未发酵杉木培养比较,DEGs主要富集在氧化还原酶活性(oxidoreductase activity);未发酵蔗渣培养与松木培养比较,DEGs主要富集在单加氧酶活性(monooxygenase activity)、氧化还原酶活性,作用于配对供体,掺入或减少氧分子(oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen)和铁离子结合(iron ion binding)等;发酵蔗渣培养与未发酵蔗渣培养比较,DEGs主要富集在单加氧酶活性、铁离子结合、次级代谢产物合成(secondary metabolite biosynthetic process)和碳水化合物代谢(carbohydrate metabolic process)等生物学过程;发酵杉木培养与松木培养比较,DEGs主要富集在β-甘露糖苷酶活性(beta-mannosidase activity);发酵蔗渣培养与松木培养比较,DEGs富集的生物学过程有醛醇分解代谢过程(alditol catabolic process)、甘油分解代谢过程(glycerol catabolic process)和细胞多糖分解代谢过程(cellular carbohydrate catabolic process),分子功能方面主要有过渡态金属离子跨膜转运活性(transition metal ion transmembrane transporter activity)和d-苏醛糖-1-脱氢酶活性(d-threo-aldose 1-dehydrogenase activity)。

图 3 差异基因GO (A)和KEGG (B)富集分析 Figure 3 GO (A) and KEGG (B) enrichment analysis of DEGs.

差异基因的KEGG通路富集分析结果见图 3B。结果显示,杉木培养与松木培养比较,DEGs未显著富集到相关KEGG通路;发酵杉木培养与未发酵杉木培养比较,DEGs显著富集的通路为星孢菌素生物合成(staurosporine biosynthesis);未发酵蔗渣培养与松木培养比较,DEGs显著富集的通路为果糖和甘露糖代谢(fructose and mannose metabolism);发酵蔗渣培养与未发酵蔗渣培养比较,DEGs显著富集的通路有戊糖和葡萄糖醛酸转换(pentose and glucuronate interconversions)、甲烷代谢(methane metabolism)和糖酵解/糖原异生(glycolysis/gluconeogenesis)等;发酵杉木培养与松木培养比较,DEGs显著富集的通路为戊糖和葡萄糖醛酸转换;发酵蔗渣培养与松木培养比较,DEGs显著富集的通路有戊糖和葡萄糖醛酸转换、半乳糖代谢(galactose metabolism)、乙醛酸盐和二羧酸盐代谢(glyoxylate and dicarboxylate metabolism)、甘油酯代谢(glycerolipid metabolism)和氰基氨基酸代谢(cyanoamino acid metabolism)等。

2.4 不同比较组差异显著基因分析

根据不同比较组间显著差异表达基因的变化倍数,分别筛选出各组间上下调排名前5的差异基因,结果如图 4A所示。从图 4A中可以看出,松木培养下表达量较高的基因主要编码MFS通用底物转运蛋白、短链脱氢酶和乙酰辅酶A合成酶类似蛋白等;杉木培养下表达量较高的基因主要编码β-半乳糖苷酶(beta-galactosidase)、脯氨酸亚氨基肽酶(proline iminopeptidase)和铁还原酶(iron reductase)等;未发酵蔗渣培养下表达量较高的基因大都功能未知;发酵蔗渣培养下表达量较高的基因主要编码己糖转运蛋白2 (hexose transporter 2)、d-半乳糖醛酸还原酶(d-galacturonate reductase)、醇酮酸还原异构酶(ketol-acid reductoisomerase)、胃蛋白酶polyporopepsin、短链脱氢酶、烯醇化酶(enolase)和SnodProt1蛋白等,而多组间共有差异基因中的锌型乙醇脱氢酶类似蛋白基因和芳基醇脱氢酶基因均在发酵蔗渣培养下表达量最低;发酵杉木培养下表达量较高的基因主要编码醛还原酶1 (aldehyde reductase 1)、甲基转移酶ausD (methyltransferase ausD)和MFS通用底物转运蛋白等。在分别以未发酵处理的松木、杉木和蔗渣培养下,CUF1-依赖的铜转运蛋白1 (CUF1-dependent copper transporter 1)、脯氨酸亚氨基肽酶、铁还原酶、黄素腺嘌呤二核苷酸-依赖的单氧酶[flavin adenine dinucleotide (FAD)-dependent monooxygenase]、乙酰辅酶A合成酶类似蛋白和细胞色素B561 (cytochrom B561)编码基因表达量均相对较高。

图 4 差异显著基因热图聚类分析 Figure 4 Heat map clustering analysis of DEGs. A: The top 5 DEGs (up and down regulated) among all six comparison groups. B: CAZymes in DEGs.

CAZymes在木质纤维素降解过程中发挥着重要作用,通过dbCAN3[20]在线注释工具对不同比较组间显著差异表达基因进行了功能注释,共注释到49个CAZymes。热图分析结果(图 4B)显示,在不同木质纤维素基质诱导培养下,广叶绣球菌CAZymes基因表达差异较大,松木培养和杉木培养表达模式较类似。亚分支1中主要为糖苷水解酶(glycoside hydrolase, GH),且在发酵甘蔗渣培养下总体表达量要高于其他基质。亚分支7中碳水化合物辅助酶(auxiliary activity, AA) AA1_1、AA3_3、AA3_2以及糖苷水解酶GH5_31和糖基转移酶(glycosyltransferase, GT) GT3编码基因在甘蔗渣培养下表达量最高。亚分支4中GH53和碳水化合物结合模块(carbohydrate-binding module, CBM) CBM21编码基因在发酵杉木培养下表达量最高。杉木培养和松木培养下GT2、AA14、AA1_2和GH18等酶的基因表达量较高。

2.5 加权基因共表达网络分析

基于不同基质诱导下的基因表达量,利用美吉生物云平台WGCNA分析功能构建加权基因共表达网络,鉴定表达模式相似的基因集合模块(module),将模块与不同木质纤维素基质类型联系起来,分析基因网络与木质纤维素基质类型之间的联系。根据计算,选择最佳软阈值(soft threshold)为7,构建无尺度共表达网络。根据分析,本研究共获得10个模块,不同模块通过不同的颜色来区分(图 5),其中turquoise模块包含的基因数目最多(1 632个),magenta模块包含的基因最少(38个)。相关性分析结果表明,green模块与未发酵蔗渣诱导培养显著正相关,相关系数为0.911 (P=0);blue模块与发酵甘蔗渣诱导培养显著正相关,相关系数为0.809 (P=0.000 26);magenta和turquoise模块与发酵杉木诱导培养显著正相关,相关系数分别为0.836 (P=0.000 1)和0.8 (P=0.000 34)。

图 5 模块与不同基质诱导相关性热图 Figure 5 The correlation heat map for modules and different substrate treatments. The left panel shows the ten modules and gene number contained in each module. The number in the rectangular box represents the correlation coefficient and corresponding P-value (in parentheses) between the module and the trait. Red block represents the positive correlation, and blue represents the negative correlation. The module with positive correlation is defined if the correlation coefficient is greater than 0.8 and P<0.05.

2.6 显著关联模块的GO和KEGG富集分析

为进一步发掘与不同基质诱导相关目的模块内基因的功能,分别对blue、green、turquoise和magenta模块内基因进行了功能富集分析。GO富集分析(图 6A)结果显示,turquoise模块内基因显著富集到尿素跨膜转运子活性(urea transmembrane transporter activity)、甲基转移酶活性(O-methyltransferase activity)和单加酶活性(monooxygenase activity)等;blue模块内基因显著富集到水解酶活性、β-甘露糖苷酶活性(beta-mannosidase activity)和3-磷酸甘油醛脱氢酶活性(glycerol-3-phosphate dehydrogenase activity)等;其余模块未检测到显著富集的GO条目。

图 6 显著相关模块基因的GO (A)和KEGG (B)功能富集 Figure 6 GO (A) and KEGG (B) enrichment analysis of genes in significantly correlated modules.

KEGG通路富集分析(图 6B)结果显示,green模块内基因显著富集的通路为星孢菌素生物合成;blue模块内基因显著富集的通路有半乳糖代谢(galactose metabolism)、果糖和甘露糖代谢(fructose and mannose metabolism)、苯丙氨酸代谢(phenylalanine metabolism)、精氨酸和脯氨酸代谢(arginine and proline metabolism)、戊糖和葡萄糖醛酸相互转化和糖鞘脂生物合成-球系列和异球系列(glycosphingolipid biosynthesis-globo and isoglobo series)。其余模块未显著富集到相关通路。

2.7 核心基因筛选

连通性表示某个基因与其他基因的关联程度,连通性越高,说明该基因在互作网络中越处于重要地位。选取blue、green、turquoise和magenta模块内连通性排序前30的基因,利用Cytoscape构建互作网络图(图 7),分析网络中的核心基因。结果显示,green模块中的EVM0003396、EVM0005606和EVM0005743,blue模块中的EVM0011996、EVM0010970和EVM0005670,turquoise模块中的EVM0001970、EVM0003293和EVM0010892,magenta模块中的EVM0000632、EVM0004451和EVM0005532分别处于各自互作网络的核心位置。核心基因功能注释信息见表 2

图 7 模块中关键基因的互作网络 Figure 7 Gene interaction network of hub genes. A: Green module. B: Blue module. C: Turquoise module. D: Magenta module.

表 2. 模块核心基因功能注释 Table 2. Functional annotation of hub genes in different modules
Module Gene ID Gene annotation
Green EVM0003396 Ring finger protein
EVM0005606 Hypothetical protein
EVM0005743 Aspyridones efflux protein
Blue EVM0011996 Hypothetical protein
EVM0010970 Predicted protein
EVM0005670 Aldo-keto reductase
Turquoise EVM0001970 Hypothetical protein
EVM0003293 Protein-lysine N-methyltransferase
EVM0010892 tRNA (His) guanylyltransferase
Magenta EVM0000632 Abhydrolase-domain-containing protein 6
EVM0004451 2-oxoglutarate-dependent dioxygenase
EVM0005532 Abhydrolase-domain-containing protein 6

3 讨论与结论 3.1 不同基质类型对基因差异表达的影响

松木屑由于松脂含量较高,大多数栽培食用菌都难以利用新鲜松木屑作为基质[12]。而广叶绣球菌栽培则可利用新鲜的松木屑作为基质,因此其对松木屑有机物质的降解相对于其他类型木质纤维素可能有独特的降解机制。杉木与松木同为针叶树软木,同样具有较高含量的挥发性分子,本研究结果也显示,松木培养与杉木培养的比较组间差异基因最少(20个,变化倍数均小于3倍),而松木培养与属于草本植物木质纤维素的蔗渣培养相比,差异基因最多。同时,相关性分析也表明基质中木质素、纤维素及半纤维素的含量影响着部分转运蛋白与脱氢酶等降解相关基因的表达。因此,广叶绣球菌对不同类型木质纤维素基质的利用采用不同的降解策略。

在以发酵蔗渣为碳源和以未发酵蔗渣为碳源培养的比较组中,单加氧酶活性显著富集,涉及的21个基因主要编码细胞色素P450单加氧酶(cytochrome P450 monoxygenase),且相对于未发酵蔗渣培养,在发酵蔗渣培养下表达量均显著下调,其中下调倍数大于3倍的包括苯甲酸4-单加氧酶(benzoate 4-monooxygenase) 编码基因(EVM0002855)和细胞色素P450单加氧酶编码基因(EVM0011015和EVM0000042)。细胞色素P450单加氧酶在木质素降解和芬顿反应中发挥重要作用[7, 21],且相对简单碳源,结构复杂的木质纤维素诱导下细胞色素P450单加氧酶均会有上调表达[7-8, 22]。发酵处理后的蔗渣由于长期堆放在自然环境下受到微生物的降解,其木质素结构受到了一定程度的破坏,因此相对于未发酵过的蔗渣细胞色素P450单加氧酶基因表达量降低。褐腐菌由于缺乏木质素降解的过氧化物酶类,其在木质纤维素解聚或修饰过程中主要依赖辅助酶类AAs[23]或芬顿反应过程产生的活性氧[24]破坏木质素结构。因此,与芬顿反应相关的铁还原酶、细胞色素B561、FAD-单加氧酶和铜转运蛋白编码基因在未发酵处理的木质纤维素样品中表达量较高。辅助酶AAs大多在未发酵处理的松木、杉木和蔗渣中表达量较高,且具有不同的诱导模式,如AA1_1、AA3_2和AA3_3等家族的酶类在未发酵蔗渣中表达量较高,而AA1_2家族的酶类在松木和杉木中表达量较高。因此,不同木质纤维素类型可能由于木质素的含量、结构等因素诱导了绣球菌不同酶系参与木质素解聚或修饰。发酵蔗渣和发酵杉木中与糖酵解相关的烯醇化酶、果胶降解相关的d-半乳糖醛酸还原酶[25]、糖转运相关的己糖转运蛋白、具有促进纤维素酶活性的胃蛋白酶polyporopepsin[26]以及通过破坏纤维素非共价键促进纤维素酶降解活性的SnodProt1蛋白[27]等蛋白的编码基因表达量较高。同时,与纤维素降解相关的GH1、GH3、GH12和GH16等家族的酶类,与半纤维素降解相关的GH2、GH27和GH31等家族的酶类以及与几丁质降解相关的GH18家族的酶类在发酵基质培养的转录组中表达量均高于未发酵处理基质培养的转录组。因此栽培基质的预处理可改善木质纤维素的理化性质,减少木质素对纤维素、半纤维素降解的结构阻碍,进而促进水解酶等降解相关基因的表达,提高降解效率。生产上可尝试通过添加发酵后的杉木或蔗渣,改善栽培菌包透气性和持水性的同时,有利于基质纤维素降解酶系的诱导表达,促进菌丝生长。木质纤维素预处理过程中会产生抑制真菌生长的醛类,因此醛还原酶在发酵杉木和发酵蔗渣为碳源的培养中表达量增高,以利于将醛类代谢为毒性更低或无毒的醇类[28]

3.2 关键模块与核心基因分析

为进一步分析广叶绣球菌利用不同基质特有的降解机制,利用WGCNA分析了基因表达模式与木质纤维素类型的相关性,通过共表达网络挖掘核心基因。与发酵蔗渣培养相关的blue模块中基因显著富集的活性有水解酶活性,这与不同基质中差异基因富集情况类似,且模块中半乳糖代谢、果糖和甘露糖代谢通路显著富集,进一步说明了纤维素和半纤维素代谢相关基因易受发酵蔗渣的诱导。而模块中的醛酮还原酶可能与木质素的解聚和修饰相关[22]。同时,甘油三磷酸脱氢酶活性在blue模块中也显著富集,甘油三磷酸脱氢酶与甘油积累相关[29],甘油调控着真菌抵抗高渗透压胁迫[30]。木质纤维素降解过程中会产生大量的乙酸盐[31-32],过量的乙酸盐能够抑制真菌的生长[33],甘油三磷酸脱氢酶基因表达量升高,可能与发酵蔗渣中水解酶基因高表达量或发酵过程中产生的乙酸盐胁迫有关。乙酸盐在乙酰辅酶A合成酶的作用下可以转化为乙酰辅酶A,过表达乙酰辅酶A合成酶可增加强微生物对乙酸盐胁迫的耐受能力[33]。不同木质纤维素为碳源培养的比较组共有的差异基因分析结果显示,乙酰辅酶A合成酶为共有的差异表达基因,推测甘油三磷酸脱氢酶和乙酰辅酶A合成酶在广叶绣球菌木质纤维素降解过程中,调控渗透压的胁迫以及木质素分解产生的乙酸盐的利用。

与未发酵蔗渣培养相关的green模块中,核心基因编码的环指蛋白(ring finger protein)属于锌指蛋白家族。根据人工锌指蛋白调控里氏木霉纤维素酶的合成[34],推测该环指蛋白在基因表达调控中起到重要作用,其在绣球菌基质降解过程中,环指蛋白可能参与了相关降解基因的表达调控。与发酵杉木培养相关的magenta模块中,核心基因编码的α/β水解酶结构域包含蛋白6 (abhydrolase-domain-containing protein 6, ABHD6)在调控脂类代谢中发挥作用,ABHD6可水解单酰基甘油和溶血磷脂酰胆碱[35],广叶绣球菌利用松木屑过程中脂肪酶(EVM0003260和EVM0000539)与其降解特性显著相关[12]。脂质代谢过程可能在广叶绣球菌木质纤维素降解过程中起到重要作用。α-酮戊二酸依赖型双加氧化酶(α-oxoglutarate-dependent dioxygenase)具有催化底物的去饱和、重排和去甲基化作用[36]。推测magenta模块中的α-酮戊二酸依赖型双加氧化酶可能与木质素的解聚和修饰相关。与发酵杉木相关的turquoise模块中,核心基因编码的tRNA (组氨酸)鸟苷转移酶[tRNA (His) guanylyltransferase]在水稻翻译效率调控和高温胁迫响应方面发挥重要作用[37],蛋白-赖氨酸N-甲基转移酶(protein-lysine N-methyltransferase)可催化赖氨酸甲基化,是一种调节多种信号通路的翻译后修饰机制[38]。tRNA (组氨酸)鸟苷转移酶和蛋白-赖氨酸N-甲基转移酶可能参与了绣球菌木质纤维素降解酶相关基因的调控。

本研究结果初步表明针对不同类型木质纤维素,广叶绣球菌采用不同的降解策略。未发酵处理的木质纤维素为基质时,与木质素结构破坏相关的基因表达量较高。经发酵处理的基质可诱导广叶绣球菌纤维素、半纤维素降解相关水解酶基因表达量升高,因此推测栽培基质的预处理对转化效率及绣球菌生长存在促进作用。WGCNA分析获得了部分与基质特异性降解相关的关键基因,相关基因可能参与了基质降解和相关基因表达的调控。同时,从研究结果也可以看到,由于广叶绣球菌基因组注释信息还不完善,部分差异显著或与降解相关的关键基因的功能还不能分析。

综合以上相关结果和分析,本研究为广叶绣球菌木质纤维素降解关键基因、分子机制研究以及栽培基质优化提供了参考。

References
[1] SISTA KAMESHWAR AK, QIN WS. Systematic metadata analysis of brown rot fungi gene expression data reveals the genes involved in Fenton's reaction and wood decay process[J]. Mycology, 2020, 11(1): 22-37 DOI:10.1080/21501203.2019.1703052.
[2] LUNDELL TK, MÄKELÄ MR, HILDÉN K. Lignin-modifying enzymes in filamentous basidiomycetes-ecological, functional and phylogenetic review[J]. Journal of Basic Microbiology, 2010, 50(1): 5-20 DOI:10.1002/jobm.200900338.
[3] HIBBETT DS, DONOGHUE MJ. Analysis of character correlations among wood decay mechanisms, mating systems, and substrate ranges in homobasidiomycetes[J]. Systematic Biology, 2001, 50(2): 215-242 DOI:10.1080/10635150151125879.
[4] KRAH FS, BÄSSLER C, HEIBL C, SOGHIGIAN J, SCHAEFER H, HIBBETT DS. Evolutionary dynamics of host specialization in wood-decay fungi[J]. BMC Evolutionary Biology, 2018, 18(1): 119 DOI:10.1186/s12862-018-1229-7.
[5] FÜCHTNER S, BROCK-NANNESTAD T, SMEDS A, FREDRIKSSON M, PILGÅRD A, THYGESEN LG. Hydrophobic and hydrophilic extractives in Norway spruce and Kurile larch and their role in brown-rot degradation[J]. Frontiers in Plant Science, 2020, 11: 855 DOI:10.3389/fpls.2020.00855.
[6] WU BJ, GASKELL J, HELD BW, TOAPANTA C, VUONG T, AHRENDT S, LIPZEN A, ZHANG JW, SCHILLING JS, MASTER E, GRIGORIEV IV, BLANCHETTE RA, CULLEN D, HIBBETT DS. Retracted and republished from "Substrate-specific differential gene expression and RNA editing in the brown rot fungus Fomitopsis pinicola"[J]. Applied and Environmental Microbiology, 2021, 87(16): e0032921 DOI:10.1128/AEM.00329-21.
[7] UMEZAWA K, NIIKURA M, KOJIMA Y, GOODELL B, YOSHIDA M. Transcriptome analysis of the brown rot fungus Gloeophyllum trabeum during lignocellulose degradation[J]. PLoS One, 2020, 15(12): e0243984 DOI:10.1371/journal.pone.0243984.
[8] GASKELL J, BLANCHETTE RA, STEWART PE, BonDURANT SS, ADAMS M, SABAT G, KERSTEN P, CULLEN D. Transcriptome and secretome analyses of the wood decay fungus Wolfiporia cocos support alternative mechanisms of lignocellulose conversion[J]. Applied and Environmental Microbiology, 2016, 82(13): 3979-3987 DOI:10.1128/AEM.00639-16.
[9] SABAT G, AHRENDT S, WU BJ, GASKELL J, HELD BW, TOAPANTA C, VUONG TV, LIPZEN A, ZHANG JW, SCHILLING JS, MASTER E, GRIGORIEV IV, BLANCHETTE RA, HIBBETT DS, BHATNAGAR J, CULLEN D. Proteome of the wood decay fungus Fomitopsis pinicola is altered by substrate[J]. Microbiology Resource Announcements, 2022, 11(9): e0058622 DOI:10.1128/mra.00586-22.
[10] LEE DJ, JANG MC, RA JO A, CHOI HJ, KIM KS, CHI YT. Noble strain of Sparassis latifolia produces high content of β-glucan[J]. Asian Pacific Journal of Tropical Biomedicine, 2015, 5(8): 629-635 DOI:10.1016/j.apjtb.2015.05.008.
[11] 肖冬来, 马璐, 杨驰, 林衍铨. 不同碳源条件下广叶绣球菌转录组分析[J]. 微生物学通报, 2019, 46(7): 1654-1661. DOI:10.13344/j.microbiol.china.180637
XIAO DL, MA L, YANG C, LIN YQ. Transcriptome analysis of Sparassis latifolia cultivated with different carbon sources[J]. Microbiology China, 2019, 46(7): 1654-1661 (in Chinese).
[12] YANG C, MA L, XIAO DL, LIU XY, JIANG XL, LIN YQ. Comparative transcriptomics reveals unique pine wood decay strategies in the Sparassis latifolia[J]. Scientific Reports, 2022, 12: 19875 DOI:10.1038/s41598-022-24171-z.
[13] 黄毅. 食用菌工厂化栽培实践[M]. 第一版. 福州: 福建科学技术出版社, 2014.
HUANG Y. Practical mushroom farming[M]. First edition. Fuzhou: Fujian Science & Technology Publishing House, 2014 (in Chinese).
[14] FOONG SY, LIEW RK, LEE CL, TAN WP, PENG WX, SONNE C, TSANG YF, LAM SS. Strategic hazard mitigation of waste furniture boards via pyrolysis: Pyrolysis behavior, mechanisms, and value-added products[J]. Journal of Hazardous Materials, 2022, 421: 126774 DOI:10.1016/j.jhazmat.2021.126774.
[15] CHEN SF, ZHOU YQ, CHEN YR, GU J. Fastp: an ultra-fast all-in-one FASTQ preprocessor[J]. Bioinformatics, 2018, 34(17): i884-i890 DOI:10.1093/bioinformatics/bty560.
[16] YANG C, MA L, XIAO DL, LIU XY, JIANG XL, YING ZH, LIN YQ. Chromosome-scale assembly of the Sparassis latifolia genome obtained using long-read and Hi-C sequencing[J]. G3:Genes Genomes Genetics, 2021, 11(8): jkab173 DOI:10.1093/g3journal/jkab173.
[17] LI B, DEWEY CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome[J]. BMC Bioinformatics, 2011, 12: 323 DOI:10.1186/1471-2105-12-323.
[18] LOVE MI, HUBER W, ANDERS S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J]. Genome Biology, 2014, 15(12): 550 DOI:10.1186/s13059-014-0550-8.
[19] XIE C, MAO XZ, HUANG JJ, DING Y, WU JM, DONG S, KONG L, GAO G, LI CY, WEI LP. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases[J]. Nucleic Acids Research, 2011, 39(suppl_2): W316-W322 DOI:10.1093/nar/gkr483.
[20] ZHANG H, YOHE T, HUANG L, ENTWISTLE S, WU PZ, YANG ZL, BUSK PK, XU Y, YIN YB. dbCAN2: a meta server for automated carbohydrate-active enzyme annotation[J]. Nucleic Acids Research, 2018, 46(W1): W95-W101 DOI:10.1093/nar/gky418.
[21] SYED K, YADAV JS. P450 monooxygenases (P450ome) of the model white rot fungus Phanerochaete chrysosporium[J]. Critical Reviews in Microbiology, 2012, 38(4): 339-363 DOI:10.3109/1040841X.2012.682050.
[22] KAMESHWAR AKS, QIN WS. Molecular networks of Postia placenta involved in degradation of lignocellulosic biomass revealed from metadata analysis of open access gene expression data[J]. International Journal of Biological Sciences, 2018, 14(3): 237-252 DOI:10.7150/ijbs.22868.
[23] RILEY R, SALAMOV AA, BROWN DW, NAGY LG, FLOUDAS D, HELD BW, LEVASSEUR A, LOMBARD V, MORIN E, OTILLAR R, LINDQUIST EA, SUN H, LaBUTTI KM, SCHMUTZ J, JABBOUR D, LUO H, BAKER SE, PISABARRO AG, WALTON JD, BLANCHETTE RA, et al. Extensive sampling of basidiomycete genomes demonstrates inadequacy of the white-rot/brown-rot paradigm for wood decay fungi[J]. Proceedings of the National Academy of Sciences of the United States of America, 2014, 111(27): 9923-9928.
[24] CASTAÑO JD, ZHANG JW, ANDERSON CE, SCHILLING JS. Oxidative damage control during decay of wood by brown rot fungus using oxygen radicals[J]. Applied and Environmental Microbiology, 2018, 84(22): e01937-18.
[25] CHROUMPI T, MÄKELÄ MR, de VRIES RP. Engineering of primary carbon metabolism in filamentous fungi[J]. Biotechnology Advances, 2020, 43: 107551 DOI:10.1016/j.biotechadv.2020.107551.
[26] SALVACHÚA D, MARTÍNEZ AT, TIEN M, LÓPEZ-LUCENDO MF, GARCÍA F, de LOS RÍOS V, MARTÍNEZ MJ, PRIETO A. Differential proteomic analysis of the secretome of Irpex lacteus and other white-rot fungi during wheat straw pretreatment[J]. Biotechnology for Biofuels, 2013, 6(1): 115 DOI:10.1186/1754-6834-6-115.
[27] PITOCCHI R, CICATIELLO P, BIROLO L, PISCITELLI A, BOVIO E, VARESE GC, GIARDINA P. Cerato-platanins from marine fungi as effective protein biosurfactants and bioemulsifiers[J]. International Journal of Molecular Sciences, 2020, 21(8): 2913 DOI:10.3390/ijms21082913.
[28] OUYANG YD, LI Q, KUANG XL, WANG HY, WU JJ, AYEPA E, CHEN H, ABRHA GT, ZHANG ZY, LI X, MA MG. YMR152W from Saccharomyces cerevisiae encoding a novel aldehyde reductase for detoxification of aldehydes derived from lignocellulosic biomass[J]. Journal of Bioscience and Bioengineering, 2021, 131(1): 39-46 DOI:10.1016/j.jbiosc.2020.09.004.
[29] ZHANG B, ZHU YL, ZHANG J, WANG DM, SUN LH, HONG J. Engineered Kluyveromyces marxianus for pyruvate production at elevated temperature with simultaneous consumption of xylose and glucose[J]. Bioresource Technology, 2017, 224: 553-562 DOI:10.1016/j.biortech.2016.11.110.
[30] SHI YK, WANG H, YAN YX, CAO HJ, LIU XH, LIN FC, LU JP. Glycerol-3-phosphate shuttle is involved in development and virulence in the rice blast fungus Pyricularia oryzae[J]. Frontiers in Plant Science, 2018, 9: 687 DOI:10.3389/fpls.2018.00687.
[31] del RÍO JC, MARQUES G, RENCORET J, MARTÍNEZ ÁT, GUTIÉRREZ A. Occurrence of naturally acetylated lignin units[J]. Journal of Agricultural and Food Chemistry, 2007, 55(14): 5461-5468 DOI:10.1021/jf0705264.
[32] WANG YN, ZHAN P, SHAO LS, ZHANG L, QING Y. Effects of inhibitors generated by dilute phosphoric acid plus steam-exploded poplar on Saccharomyces cerevisiae growth[J]. Microorganisms, 2022, 10(7): 1456 DOI:10.3390/microorganisms10071456.
[33] DING J, HOLZWARTH G, PENNER MH, PATTON-VOGT J, BAKALINSKY AT. Overexpression of acetyl-CoA synthetase in Saccharomyces cerevisiae increases acetic acid tolerance[J]. FEMS Microbiology Letters, 2015, 362(3): 1-7.
[34] 孟庆山, 李嘉祥, 张飞, 赵心清, 白凤武. 人工锌指蛋白介导调控的里氏木霉纤维素酶生产[J]. 生物工程学报, 2019, 35(1): 81-90.
MENG QS, LI JX, ZHANG F, ZHAO XQ, BAI FW. Artificial zinc finger protein mediated cellulase production in Trichoderma reesei Rut-C30[J]. Chinese Journal of Biotechnology, 2019, 35(1): 81-90 (in Chinese).
[35] LI KS, EGELANDSDAL B, OLSEN RE. Hydrolysis activity of pyloric cecal enterocytes of brown trout (Salmo trutta) toward monoacylglycerol and lysophosphatidylcholine[J]. Lipids, 2018, 53(6): 615-625 DOI:10.1002/lipd.12056.
[36] ABE I. Nonheme iron- and 2-oxoglutarate-dependent dioxygenases in fungal meroterpenoid biosynthesis[J]. Chemical and Pharmaceutical Bulletin, 2020, 68(9): 823-831 DOI:10.1248/cpb.c20-00360.
[37] CHEN K, GUO T, LI XM, ZHANG YM, YANG YB, YE WW, DONG NQ, SHI CL, KAN Y, XIANG YH, ZHANG H, LI YC, GAO JP, HUANG XH, ZHAO Q, HAN B, SHAN JX, LIN HX. Translational regulation of plant response to high temperature by a dual-function tRNAHis guanylyltransferase in rice[J]. Molecular Plant, 2019, 12(8): 1123-1142 DOI:10.1016/j.molp.2019.04.012.
[38] LEVY D. Lysine methylation signaling of non-histone proteins in the nucleus[J]. Cellular and Molecular Life Sciences, 2019, 76(15): 2873-2883 DOI:10.1007/s00018-019-03142-0.