崔 筱，张玉亭，孔维丽，胡素娟，刘 芹，王彦坡，吴 杰，师子文

（1. 河南省农业科学院 植物营养与资源环境研究所，河南 郑州 450002；2. 西峡县食用菌生产办公室，河南 西峡 474599）

平菇（Pleurotus ostreatus）是一种常见的药食两用的真菌，具有抗肿瘤、抗氧化、免疫调节等多种药理活性，已成为我国栽培最为广泛的食用菌之一[1-2]。平菇菌丝的生长容易受到生长环境，尤其是外源调节剂的影响，进而影响其子实体发育过程[3-4]。生长素（Auxin）是一种植物激素，吲哚乙酸（IAA）是一种重要的生长素，是细胞增殖与伸长的重要调节剂，对细胞的周期进展及相关蛋白质也具有重要的调控作用。研究表明，外源IAA 能够反馈调节平菇内源IAA 的分泌及生长素氧化酶的活性，影响菌丝的生长速度及平菇的生长发育过程[5]，且通过转录组分析证明，几丁质代谢、色氨酸代谢和MAPK 信号通路是IAA 影响平菇菌丝生长的潜在分子机制[6]。而子实体作为平菇主要的食用部分，外源IAA 对其生长发育的影响机制研究未见报道，但转录组学技术在食用菌的子实体发育及环境响应机制研究中已被广泛应用。研究子实体对外源IAA响应的机制，对平菇子实体阶段的生长管理具有重要意义。鉴于此，设置高、低浓度外源IAA 处理不同生育阶段的平菇，通过转录组技术筛选平菇不同生育阶段对IAA 响应的差异表达基因（DEGs），建立基因共表达网络，进而揭示IAA 对平菇生长发育的影响机制，为外源IAA的施用提供理论指导。

1 材料和方法

1.1 材料

供试菌株：平菇P99，由河南省食用菌种质资源库提供。供试母种培养基：马铃薯葡萄糖琼脂培养基（PDA），购自北京奥博星生物技术有限责任公司。培养料配方：棉籽壳93%、麸皮5%、石灰2%。IAA购自Solarbio公司。

1.2 方法

1.2.1 不同浓度IAA 的配制 利用无水乙醇作为溶剂，配制成浓度分别为1×10-3（T3）、1×10-8（T8）、0（CK）mol/L 的IAA 溶液，用0.22 μm 微孔滤膜过滤除菌后储存于4 ℃冰箱备用。

1.2.2 样品制备 菌丝体收集按照崔筱等[6]的方法。用直径为5 mm 的打孔器将活化后的平菇菌株分别接种到直径为15 cm、含有3 种不同浓度（T3、T8、CK）IAA 溶液的PDA 平板上，置于25 ℃恒温培养箱中黑暗培养7 d或15 d。试验设置6组处理，分别记为T37d、T315d、T87d、T815d、CK7d 和CK15d，每处理重复3 次，每个重复3 皿。取不同处理的菌丝体0.5 g，分别装入5 mL 的无菌冻存管中，液氮中冷冻，于-80 ℃冰箱保存备用。子实体收集方法具体为：选用14 cm×28 cm×0.005 cm聚丙烯栽培袋，每袋质量250 kg，干料质量150 kg，菌丝满袋后，开口出菇，每袋喷施10 mL 不同浓度（T3、T8）的IAA 溶液，以不喷施生长素的菌袋作为CK，设置3个重复。喷施IAA 后，保持菇房内温度15～20 ℃，相对湿度80%～90%，二氧化碳浓度1 178.57～1 571.42 mg/m3，光照强度150 lx，取不同处理的25 d、27 d、31 d 的样品0.5 g，分别装入5 mL的无菌冻存管中，液氮中冷冻，于-80 ℃冰箱保存备用。

1.2.3 总RNA 的提取、文库构建和转录组测序 按照崔筱等[6]的方法，利用Trizol reagent（Takara）法提取各处理样品总RNA，并将提取的不同重复RNA进行等量混合，Nanodrop 2000分光光度计（美国，赛默飞）对所提RNA 的浓度和纯度进行检测，琼脂糖凝胶电泳检测RNA 完整性，合格后送样至上海美吉生物医药科技有限公司进行转录组测序。

1.2.4 数据获取 研究数据为平菇组织的转录组原始测序数据（https://cloud.majorbio.com/）。共包含了15 个样本，分别是CK 组（CK7d、CK15d、CK25d、CK27d、CK31d）、高浓度T3 处理组（T37d、T315d、T325d、T327d、T331d）、低浓度T8 处理组（T87d、T815d、T825d、T827d、T831d）。其中7 d、15 d 是菌丝体阶段，25 d、27 d、31 d 是子实体阶段，25 d 是子实体发育阶段中的原基分化期。

1.2.5 DEGs 筛选 利用DESeq2[7]（Version 1.34.0，https://bioconductor.org/packages/release/bioc/html/DESeq2.html）按照时间梯度分别对CK7d vs T37d、CK15d vs T315d、CK25d vs T325d、CK27d vs T327d、CK31d vs T331d 及CK7d vs T87d、CK15d vs T815d、CK25d vs T825d、CK27d vs T827d、CK31d vs T831d 的基因进行差异分析，所有基因经过分析得到相应的P值和lgFC值，阈值设定：P值＜0.05，且|logFC|＞1。

通过aov 函数使用F检验方法对相同浓度（T3、CK、T8）下不同时间点样本（7 d、15 d、25 d、27 d、31 d）间的基因进行差异分析，得到相应的P值，采用Benjamini&Hochberg 方法进行多重检验校正，得到校正的P值即Q值，DEGs阈值设定：Q值＜0.01。

1.2.6 Venn 分析 利用在线工具Venn[8]（http://bioinformatics.psb.ugent.be/webtools/Venn/）分 别 对CK7d vs T37d 的DEGs 与CK15d vs T315d 的DEGs、CK7d vs T87d 的DEGs 与CK15d vs T815d 的DEGs取交集，分别获得菌丝体的高浓度和低浓度下的交集DEGs。

同样利用在线工具Venn 分别对CK25d vs T325d 的DEGs、CK27d vs T327d 的DEGs 与CK31d vs T331d 的DEGs、CK25d vs T825d 的DEGs、CK27d vs T827d 的DEGs 与CK31d vs T831d 的DEGs 取 交集，分别获得子实体的高浓度和低浓度下的交集DEGs。

1.2.7 DEGs 表达趋势聚类分析 基于上述时间梯度得到10 组DEGs，对相同时间不同IAA 浓度下的DEGs 取并集，共得到5 组（7 d、15 d、25 d、27 d、31 d）并集DEGs，结合各自相同时间下RNA-seq 数据表达矩阵，得到DEGs 在不同浓度（CK、T8、T3）下的表达矩阵。分别将5 组DEGs 表达矩阵输入到Stem 软件[9]中进行表达趋势聚类分析，通过参数设置（Clustering method：STEM Clustering Method；Significance level：0.05），分别进行5 次聚类分析后，得到相同时间不同浓度下DEGs 的不同模块基因集。

基于上述相同浓度不同时间点样本得到的3组DEGs，并结合各自相同浓度下RNA-seq 数据表达矩阵，得到DEGs在不同时间点（7 d、15 d、25 d、27 d、31 d）下的表达矩阵。分别将3 组DEGs 表达矩阵输入到Stem 软件中进行表达趋势聚类分析，通过参数设 置（Clustering method：STEM Clustering Method；Significance level：0.05），分别进行3次聚类分析后，得到相同浓度不同时间点样本间DEGs 的不同模块基因集。

1.2.8 DEGs 功能富集分析 基因本体论（Gene ontology，GO）系统包括3 个部分：生物学过程（Biological process，BP）、分 子 功 能（Molecular functions，MF）、细胞组分（Cellular components，CC）。利 用 平 台 软 件Goatools[10]（https：//github. com/tanghaibao/GOatools）分别对菌丝体下高浓度间DEGs 和低浓度间DEGs、子实体下高浓度间DEGs和低浓度间DEGs 参与的GO 系统和KEGG 通路进行富集分析。显著性阈值P值＜0.05 视为显著富集结果，选择排名前10（按照P值排名，若富集结果不足10 条则进行全部展示）的GO 注释和KEGG 通路注释进行展示。

1.2.9 DEGs 共表达分析 基于上述得到的菌丝体高浓度DEGs、菌丝体低浓度DEGs、子实体高浓度DEGs 以及子实体低浓度DEGs，结合基因的表达量，分别计算出菌丝体高浓度DEGs、菌丝体低浓度DEGs、子实体高浓度DEGs 以及子实体低浓度DEGs 间的Pearson 相关性系数，得到相应的P值和相关性系数cor，并进行相关性检验。此外，采用“BH”校正，得到校正后的Q值，根据阈值Q值＜0.05，且|cor|＞0.5，得到基因共表达关系对。接着利用Cytoscape 软件[11]进行共表达网络构建（若共表达关系对过多，按照相关性系数cor 排名，分别选取正负相关各排名前50的关系对进行网络构建）。

2 结果与分析

2.1 DEGs统计结果

根据平菇生育阶段的变化，按照所设阈值筛选低浓度与高浓度外源IAA 条件下不同时间点的DEGs（图1），具体数目如表1所示。

表1 不同生育阶段平菇对外源IAA浓度响应的DEGs数量统计Tab.1 The number of DEGs at different growth stages of P.ostreatus in response to the concentration of exogenous IAA

图1 T3（A）和T8（B）IAA浓度下不同时间点的DEGs火山图Fig.1 Volcano plots of DEGs at different time points under T3（A）and T8（B）IAA concentrations

菌丝体阶段（7 d 与15 d），T3 浓度下共筛选到DEGs 2 931 个，其中上调DEGs 1 593 个，下调DEGs 1 338个；而T8浓度下共筛选到DEGs 2 007个，包括上调DEGs 713 个，下调DEGs 1 294 个。子实体阶段，T3 浓度下共筛选到DEGs 3 978 个，其中上调DEGs 2 284 个，下调DEGs 1 694 个；T8 浓度下DEGs共3 956 个，包括上调DEGs 1 894 个，下调DEGs 2 062 个。此外，在相同浓度下的不同生育阶段进行DEGs 分析，结果显示，CK、T8、T3 浓度下不同时间点样本间DEGs分别为3 209、2 856、3 173个。

2.2 Venn分析结果

如图2A 所示，T3 浓度下，菌丝体阶段的交集DEGs 有217 个。其中包括上调表达的TRINITY_DN8187_c0_g1、TRINITY_DN572_c3_g3、TRINITY_DN1371_c0_g1等及下调表达的TRINITY_DN14069_c0_g1、TRINITY_DN14570_c0_g2、TRINITY_DN41591_c0_g1、TRINITY_DN27347_c0_g1、TRINITY_DN35150_c0_g1 等。而T8 浓度下，菌丝体阶段的交集DEGs为81 个，结果如图2B 所示。其中包括上调表达的TRINITY_DN2732_c0_g1、TRINITY_DN8088_c0_g1、TRINITY_DN9337_c0_g1 等及下调表达的TRINITY_DN5232_c1_g1、TRINITY_DB2618_c0_g1、TRINITY_DN601_c0_g1等。

图2 DEGs Venn图Fig.2 Venn diagram of DEGs

子实体阶段，T3与T8浓度的交集DEGs分别为10个与87个，如图2C—D所示。其中，T3浓度下包括上调 表 达 的 TRINITY_DN4505_c0_g1、TRINITY_DN7383_c0_g1 等及下调表达的TRINITY_DN25720_c0_g1、TRINITY_DN18000_c0_g1、TRINITY_DN13805_c0_g1、TRINITY_DN5144_c0_g1等。T8浓度下包括上调 表 达 的 TRINITY_DN1089_c0_g2、TRINITY_DN181_c0_g1、TRINITY_DN3832_c0_g1等及下调表达的TRINITY_DN560_c0_g1、TRINITY_DN530_c0_g1、TRINITY_DN4090_c0_g1等。

2.3 DEGs表达趋势聚类分析结果

相同时间点的不同IAA 浓度下，时间点7 d 下并集DEGs 有2 358 个，聚类为16 个模块基因集（图3A）；在时间点15 d下并集DEGs有1 362个，聚类分析得到16 个模块基因集（图3B）；时间点25 d 下并集DEGs 有2 294 个，聚类分析为16 个模块基因集（图3C）；时间点27 d 下并集DEGs 有3 234 个，聚类分析得到16 个模块基因集（图3D）；时间点31 d 下并集DEGs 有905 个，聚类分析后，得到16 个模块基因集（图3E）。且5 个不同的时间点下均有符合先上升后下降或先下降后上升的趋势模块基因集，均得到2个趋势模块基因集（模块9 和模块6，其中，模块9 为先上升后下降，模块6 为先下降后上升）。其中，菌丝体交集DEGs TRINITY_DN8187_c0_g1、TRINITY_DN27347_c0_g1、TRINITY_DN601_c0_g1 同属于模块6，而TRINITY_DN14069_c0_g1归于模块9。子实体交集DEGs TRINITY_DN5144_c0_g1 及TRINITY_DN4090_c0_g1 同 属 于 模 块6，TRINITY_DN4505_c0_g1、TRINITY_DN25720_c0_g1、TRINITY_DN18000_c0_g1、TRINITY_DN13805_c0_g1 同 属 于模块9。

图3 相同时间点不同IAA浓度下Stem聚类趋势基因集Fig.3 Stem clustering trend gene set of different concentration of IAA at the same time points

如图4A—C 所示，相同IAA 浓度不同时间点下，CK、T8、T3浓度下均得到50个模块基因集，其中符合一直上升或者一直下降趋势的模块均为3 个（模块41、模块38和模块9）。

图4 相同IAA浓度不同时间点下Stem聚类趋势基因集Fig.4 Stem clustering trend gene set of the same concentraion of IAA at different time points

各处理趋势模块数量如表2 所示。其中，菌丝体交集DEGs TRINITY_DN572_c3_g3、TRINITY_DN1371_c0_g1、TRINITY_DN2732_c0_g1 同属于上升趋势模块，TRINITY_DN41591_c0_g1、TRINITY_DN35150_c0_g1、TRINITY_DN5314_c0_g2、TRINITY_DN43009_c0_g1、TRINITY_DN5232_c1_g1、TRINITY_DN2618_c0_g1 同属于下降趋势模块。子实体交集DEGs TRINITY_DN1089_c0_g1、TRINITY_DN181_c0_g1、TRINITY_DN3832_c0_g1 同属于上升趋势模块，TRINITY_DN560_c0_g1、TRINITY_DN530_c0_g1 及TRINITY_DN4090_c0_g1同属于下降趋势模块。

表2 不同处理中各趋势模块的DEGs数量统计Tab.2 The number of DEGs in each trend module with different treatments

2.4 DEGs功能富集分析结果

GO 系统分析表明，菌丝体下高浓度间DEGs 共富集得到52 个GO BP（主要涉及到溶酶体、细胞分布、细胞凋亡等生物学功能）、1 个GO CC（胞外区域）、17 个GO MF（主要为氧化还原酶活性、几丁质结合等分子功能）。KEGG 富集分析结果表明，可富集到1 条KEGG 通路（甘油磷脂代谢）（图5A）。菌丝体低浓度DEGs共富集得到58个GO BP（主要涉及到多种化合物的代谢、生物合成及分解过程）、12个GO CC（质膜、核糖体等细胞组分）、38 个GO MF（多种化合物裂解酶及解氨酶活性）和4 条KEGG 通路（苯丙氨酸、色氨酸、蛋氨酸等的代谢）（图5B）。

子实体高浓度DEGs 共富集得到10 个GO BP（主要涉及到ADP 代谢与ATP 合成以及丙酮酸代谢等生物学功能）、1 个GO CC（磷酸丙酮酸水合酶复合物）、4 个GO MF（氧化还原酶活性、FMN 结合、NAD（P）H 脱氢酶及磷酸丙酮酸水合酶活性）和4条KEGG 通路（糖酵解、RNA 降解、甲烷代谢及醌类的生物合成）（图5C）。子实体低浓度DEGs 共富集得到22 个GO BP（呼吸链的组装、膜脂分解代谢及乙醇代谢等生物学功能）、16个GO MF（过渡金属离子结合、血红素结合、氧化还原酶活性、铁结合等分子功能）和6 条KEGG 通路（过氧化物酶体、半乳糖代谢、鞘脂代谢等）（图5D）。这里按照P值排名，选取GO和KEGG 排名前10的注释进行展示。

图5 不同IAA浓度下DEGs的GO和KEGG富集结果Fig.5 GO and KEGG enrichment results of DEGs at different IAA concentrations

2.5 DEGs共表达分析结果

分别计算菌丝体高浓度DEGs、菌丝体低浓度DEGs、子实体高浓度DEGs 以及子实体低浓度DEGs 间的Pearson 相关性系数，根据阈值Q值＜0.05且|cor|＞0.5，其中菌丝体高浓度DEGs 共得到6 735个共表达关系对，菌丝体低浓度DEGs共得到766个共表达关系对，子实体高浓度DEGs 共得到32个共表达关系对，子实体低浓度DEGs 共得到1 565 个共表达关系对，分别选取正负相关排名前50的关系对进行网络构建，如图6所示。

结果表明，T3 菌丝体阶段的DEGs TRINITY_DN8187_c0_g1与多种基因表现出共表达关系，同时，TRINITY_DN572_c3_g3、 TRINITY_DN1371_c1_g1、TRINITY_DN14069_c0_g1、TRINITY_DN14570_c0_g2、TRINITY_DN41591_c0_g1、TRINITY_DN27347_c0_g1、TRINITY_DN5314_c0_g2 等 交 集DEGs 均 表现较多的共表达关系，且TRINITY_DN14570_c0_g2、TRINITY_DN41591_c0_g1、TRINITY_DN27347_c0_g1构成同一共表达网络（图6A）。

T8菌丝体阶段DEGs TRINITY_DN9079_c0_g1、TRINITY_DN5232_c1_g1、TRINITY_DN2618_c0_g1及TRINITY_DN601_c0_g1 表现出较多共表达关系，且属同一共表达网络（图6B）。

T3 子实体阶段DEGs TRINITY_DN4505_c0_g1、TRINITY_DN25720_c0_g1、TRINITY_DN18000_c0_g1、TRINITY_DN13805_c0_g1 及TRINITY_DN5144_c0_g1均表现出较多共表达关系（图6C）。

T8子实体阶段DEGs TRINITY_DN181_c0_g1、TRINITY_DN3832_c0_g1、TRINITY_DN560_c0_g1、TRINITY_DN530_c0_g1 及TRINITY_DN4090_c0_g1表现出较多共表达关系，TRINITY_DN1089_c0_g1与TRINITY_DN560_c0_g1 及TRINITY_DN530_c0_g1属同一个共表达网络（图6D）。

图6 不同浓度下DEGs共表达网络图Fig.6 The co-expression network of DEGs under different IAA concentrations

3 结论与讨论

本研究通过转录组分析平菇不同生长阶段对不同浓度外源IAA 响应的DEGs，结果显示，平菇在不同生长阶段对外源IAA 浓度的响应不同。菌丝体阶段，高浓度外源IAA 引起DEGs 数量高于低浓度外源IAA，说明菌丝体阶段平菇对高浓度外源IAA 响应更为强烈；子实体阶段，不同浓度外源IAA引起的DEGs 数量相近。菌丝体阶段，不同浓度IAA 下的交集DEGs 较多，而子实体阶段，不同浓度IAA 下的交集DEGs 较少。因此，平菇不同生长阶段对不同浓度外源IAA 的响应机制不同，下面对DEGs的功能进行进一步分析。

菌丝体是平菇的营养器官，菌丝体通过吸收利用培养基质中的纤维素、氮、磷、钾等养分，完成该阶段发育。本试验中，低浓度外源IAA 能够诱导核糖体及质膜相关基因的异常表达，以及氨基酸（如苯丙氨酸、色氨酸、蛋氨酸等）代谢相关的基因也表现出对低浓度IAA的响应。之前研究发现，外源IAA浓度的降低能够促进菌丝体的生长，低至1×10-10mol/L时，菌丝生长速度变慢[5]。核糖体是蛋白质翻译的重要场所，而苯丙氨酸、色氨酸及蛋氨酸能够影响植物激素例如乙烯、吲哚乙酸及木质素、花青素等的合成，对菌丝的生长具有重要作用[12]。之前有报道证明，2×10-4mol/L 外源IAA 能够影响桃果实与种子中生长素、茉莉酸及脱落酸含量，与内源植物激素共调控果实生长发育及成熟过程[13]。本研究发现，上调表达的TRINITY_DN2732_c0_g1 富集于核糖体GO 结果中，因此，预测在菌丝体阶段，低浓度外源IAA 通过TRINITY_DN2732_c0_g1 调控核糖体合成内源植物激素，调控平菇菌丝体阶段的生长发育。高浓度IAA 下的DEGs 主要富集于植物的氧化应激相关生物学过程，涉及到溶酶体、细胞凋亡等生物学功能，氧化还原酶、几丁质结合等分子功能及甘油磷脂代谢过程。甘油磷脂代谢的产物能够破坏细胞膜，促进细胞凋亡，衰老的细胞器或生物大分子被溶酶体吞噬，促进细胞的更新[14]。甘油磷脂代谢过程富集到下调的DEGs TRINITY_DN41591_c0_g1。几丁质作为细胞壁的主要成分，其代谢主要通过自溶过程实现，在几丁质代谢过程的富集结果中，包含了菌丝体阶段下调表达的TRINITY_DN14069_c0_g1。此外，TRINITY_DN35150_c0_g1、TRINITY_DN5314_c0_g2 及TRINITY_DN43009_c0_g1 均富集在氧化还原酶活性的分子功能，且三者均富集于下降趋势的模块，此外，三者共同构成的共表达网络结果显示，三者表达水平存在正相关性。之前针对高浓度IAA 下平菇菌丝体关键基因的挖掘也同样发现，高浓度IAA 能够引起氧化还原酶相关基因的异常表达，从而保护平菇细胞免受氧化损伤[5]。因此，高浓度的IAA 对平菇具有一定的胁迫作用，同时会抑制TRINITY_DN41591_c0_g1调控的甘油磷脂代谢过程，TRINITY_DN14069_c0_g1调控的几丁质代谢及TRINITY_DN35150_c0_g1、TRINITY_DN5314_c0_g2 及TRINITY_DN43009_c0_g1 调节的氧化还原酶活性，以提高平菇对高浓度IAA的抗性。

子实体阶段涉及到与菌丝体阶段截然不同的分子机制，低浓度IAA 下的平菇子实体阶段DEGs主要富集在代谢过程和呼吸链组分中，与之前报道平菇子实体发育阶段的基因转录组分析结果一致[15]。富集结果显示，下调表达的TRINITY_DN530_c0_g1 能够调控呼吸链的组分，上调表达的TRINITY_DN3832_c0_g1 则对鞘脂的代谢具重要作用。鞘脂作为生物膜的重要组成成分，其代谢过程及相应代谢产物对细胞的生长、分化及衰老具有重要意义，同时也可作为这些代谢过程的重要信号分子[16-17]。 因 此，子 实 体 阶 段 低 浓 度IAA 通 过TRINITY_DN530_c0_g1 及TRINITY_DN3832_c0_g1调控呼吸链组分及鞘脂代谢。高浓度IAA 引起子实体阶段ATP 合成相关基因的异常表达，例如糖酵解途径、ADP 代谢、ATP 合成等。子实体阶段平菇生长迅速，代谢速率较快，因此，伴随大量ATP 的合成及相关过程的激活。ATP 的合成除能为平菇子实体的生长提供能量外，也能够提高其对环境胁迫的抗性。例如，FENG 等[18]研究发现，真菌对铜胁迫的耐受程度与ATP 合成、糖酵解途径活性密切相关。高浓度IAA下GO及KEGG富集结果显示，上调表达的TRINITY_DN21430_c0_g1能够调控糖酵解、RNA 降解及甲烷代谢过程，进而促进平菇子实体阶段对高浓度IAA的响应。

综上所述，不同浓度IAA 对平菇菌丝体与子实体阶段的影响机制不同。低浓度IAA 能够通过TRINITY_DN2732_c0_g1 调控核糖体合成内源植物激素，同时通过 TRINITY_DN530_c0_g1 及TRINITY_DN3832_c0_g1 促进子实体阶段的代谢进程。高浓度IAA 能够通过TRINITY_DN41591_c0_g1调控甘油磷脂代谢过程，TRINITY_DN14069_c0_g1调控几丁质代谢及TRINITY_DN35150_c0_g1、TRINITY_DN5314_c0_g1、TRINITY_DN43009_c0_g1调节氧化还原酶活性，同时通过TRINITY_DN21430_c0_g1调控子实体阶段ATP的合成。