低温诱导稻曲病菌菌核形成的转录组学分析

吕楚阳，邓平川，张晓丽，孙钰超，梁五生，胡东维

（1浙江大学生物技术研究所/水稻生物学国家重点实验室，杭州 310058；2浙江大学作物科学研究所，杭州 310058）

0 引言

【研究意义】近年来，我国稻曲病已从水稻次要病害上升为主要病害，成为全世界稻米主产区的主要真菌病害之一[1]。稻曲病是由子囊菌门麦角菌科的稻曲病菌（Villosiclava virens）侵染水稻颖花造成的。病原菌侵染后，整个颖花形成球状菌落，最外层有大量黄色和墨绿色的厚垣孢子层。在一定环境条件下，稻曲球表面可形成一个至数个菌核[2]。菌核及其有性生殖是稻曲病菌安全越冬和提供来年初侵染源的重要方式，菌核及其萌发过程的控制是稻曲病控制的中心环节和瓶颈问题。因此，深入研究菌核的分化与形成机制对探索稻曲病防控新策略具有重要意义。【前人研究进展】在稻曲病的侵染循环过程中，病菌菌核常被认为是造成病害初次侵染的重要来源。稻曲病菌菌核可以安全越冬，并在来年萌发产生大量的子囊孢子，且在时间上与稻曲病菌在水稻孕穗期侵染基本一致[3-6]。稻曲病菌菌核一般形成在高海拔或温带地区，且出现在稻曲球发育的后期[7-8]，而在我国南方地区数量很少，在部分地区如广东尚未发现菌核[9]。但也有研究发现，位于长江中下游地区、属于亚热带地区的浙江省，稻曲病菌不但可产生菌核，且菌核的数量庞大[10]。笔者课题组近期研究表明，在稻曲球发育前期，低温是诱导稻曲病菌产生菌核的重要环境因子之一，室内夜间15℃处理3 d即足以诱导菌核形成[11]；减少田间菌核数量可显著减轻稻曲病的发生[12]。为了验证低温是菌核产生的主导环境因子，笔者对低温诱导后的稻曲球进行了系统性解剖，结果发现稻曲球内部具有比以往的认识更多的隐含菌核。目前，转录组测序技术已经成功运用到病原菌形态发育研究中，YU等[13]利用转录组技术比较分析稻曲病菌子实体发育和产分生孢子阶段的基因表达差异，揭示了稻曲病菌的有性生殖相关候选基因；韩彦卿等[14]采用测序技术对接种稻曲病菌的抗、感病水稻进行了转录组学分析，初步探明参与水稻与稻曲病菌之间早期互作的调控网络。【本研究切入点】通过转录组测序技术，全面、快速地获取低温诱导菌核形成过程中的所有转录本信息，并系统分析相关的基因表达变化。【拟解决的关键问题】采用转录组测序技术对低温处理的稻曲病菌样品进行研究，探明稻曲病菌菌核形成初期的基因表达变化，阐释低温诱导稻曲病菌菌核形成的分子机制并了解其发育进程，为稻曲病防控提供理论依据。

1 材料与方法

试验于2018—2019年在浙江大学完成。

图1 晚秋田间形成的稻曲球Fig. 1 The rice false smut balls formed in late autumn in paddy field

1.1 试验材料

对照组稻曲球样本采自浙江省杭州市桐庐县的稻曲病常发生地块，采集地区前10 d的日间平均气温为31.2℃。另采集浙江省宁波市奉化区前10 d的日间平均气温为23.2℃的稻曲球为处理组样本，本组样本的田间温度与对照组相差约10℃，视为低温处理（图1）。采集发育前期的稻曲球立即放入液氮中速冻，-80℃冰箱中保存。测序所用样品均设置3个生物学重复，每个重复3—5个稻曲球，分别记为对照组：TL911_1、TL911_2、TL911_3；低温处理组：FH1016_1、FH1016_2、FH1016_3。

1.2 试验方法

1.2.1 建库测序 将提取的稻曲病菌RNA利用琼脂糖凝胶电泳、Nanodrop、Qubit 2.0 Fluorometer、Agilent 2100等仪器进行降解程度、是否污染、纯度、浓度和完整性共5个方面的评估，质检合格后，以片段化的mRNA为模板，随机寡核苷酸为引物，在 M-MuLV逆转录酶体系中合成cDNA第一条链，随后用RNaseH降解 RNA链，并在 DNA polymerase I体系下，以dNTPs为原料合成 cDNA第二条链。纯化后的双链cDNA经过末端修复、加 A尾并连接测序接头，用AMPure XP beads筛选250—300 bp左右的cDNA，进行PCR扩增并再次使用AMPure XP beads纯化PCR产物，构建cDNA文库，并委托北京诺禾致源生物信息科技有限公司利用Illumina HiSeqTM4000测序平台进行高通量测序。

1.2.2 转录组数据分析 利用 Hisat2（v2.1.0）[15]将各样本数据（clean reads）分别与稻曲病菌“UV-8b”参考基因组序列（https://www.ncbi.nlm.nih.gov/genome/31935?genome_assembly_id=59177）进行比对，保留唯一基因组比对序列，结合基因组注释文件，利用StringTie（v2.0）软件完成各基因定量，差异表达基因（differentially expressed gene，DEG）鉴定依赖于R包ballgown，筛选标准为|log2fold change (FC)|≥1且q-value≤0.05[16]。采用 Gene Ontology（GO）数据库（http://geneontology.org）和 Kyoto Encyclopedia of Genes and Genomes（KEGG）数据库（https://www.genome.jp/kegg/pathway.html）对DEG分别进行GO和Pathway富集分析。

1.2.3 差异表达基因的 qRT-PCR验证 为验证转录组数据准确性，随机选取6个DEG进行qRT-PCR定量。利用Primer Premier 5.0设计qRT-PCR引物，内参基因选用α-tubulin，引物列表详见表 1。随后采用仪器CFX96 Real-Time System（BIO-RAD）进行qPCR扩增，反应程序：95℃ 5 min；95℃ 10 s，60℃ 30 s，40个循环，每个样品设置3次重复。采用2-ΔΔCT法计算差异基因相对表达量[17]。

表1 实时荧光定量PCR所用引物Table 1 Primers used in qRT-PCR

2 结果

2.1 测序数据质量分析

转录组分析共计6个样本，包括稻曲病菌对照组样本 TL911_1、TL911_2、TL911_3和低温处理组样本FH1016_1、FH1016_2、FH1016_3。过滤带有接头（adapter）序列和低质量reads，共得到59.78 G的clean data，样本平均数据量约为9.96 G，其中，近93.17%的序列能够比对到参考基因组（表 2）。以上结果表明测序产出数据质量良好，可进行后续的生物信息学分析。

基于各基因表达量，计算样本间相关系数（R2），R2越接近 1，表明样本之间表达模式的相似度越高。由图2-A可知，对照组、低温处理组3个生物学重复间任意两个样本R2均＞0.99，表明生物学重复性好、数据可靠。此外，比较了不同处理下的基因表达相关性，发现对照组TL911与处理组FH1016之间的平均R2为0.93，低于组内生物学重复，说明低温处理组的基因表达与对照组之间存在一定差异，暗示低温诱导稻曲病菌菌核形成过程中确实发生了基因转录变化。样本间层级聚类也有着相近的结果，即相同样本、不同生物学重复优先聚类（图2-B）。

表2 转录组测序数据以及与参考基因组比对结果Table 2 Sequencing data and alignment results of sequencing data with the reference genome

图2 样本间关系的整体特征Fig. 2 Overall characterization of sample relationship

综合TL911和FH1016样本基因表达分析，共检测到8 426个基因发生一定程度表达（图3-A），占稻曲病菌基因组基因的97.13%。按照DEG筛选标准，只有793个基因表达量发生了显著性变化，占总体基因的9.41%，其中，398个基因表达上调，395个基因表达下调（图 3-B）。上、下调表达倍数最高的前 6位基因及其相对表达量、FPKM表达量分别见图3-C、3-D。

2.2 GO富集分析

GO注释主要从3个角度描述基因的生物学功能，包括分子功能、生物过程和细胞成分。DEG在分子功能分类中，显著富集于氧化还原酶活性（GO:0016491），其上调基因 142个，下调基因89个；阳离子结合（GO:0043169），其上调基因146个，下调基因83个；金属离子结合（GO:0046872），其上调基因145个，下调基因83个；过渡金属离子结合（GO:0046914），其上调基因113个，下调基因 61个（图 4-A）；而在生物过程类别中，碳水化合物代谢过程（GO:0005975）和氧化还原进程（GO:0055114）等通路富集程度最高（图4-B）；进一步对分子功能中氧化还原酶活性功能条目的DEG分析，发现NADPH氧化酶（KDB12731.1）、过氧化氢酶（KDB11363.1）、FAD依赖性氧化还原酶（KDB11634.1）、肌氨酸氧化酶（KDB12820.1）在低温处理稻曲病菌中显著上调，结合生物过程类别中氧化还原过程的相对活跃，推测这些真菌机体的氧化还原反应可能与稻曲病菌菌核形成有关。另外，在所有的 GO分类中包含上调与下调基因的表达，但整体上并未发现有显著富集的下调功能条目。

图3 低温处理下稻曲病菌表达差异显著基因的比较分析Fig. 3 The comparative analysis of DEGs in V. virens under low temperature

图4 稻曲病菌低温诱导菌核形成相关基因的鉴定及富集分析Fig. 4 Identification and enrichment analysis of genes involved in low temperature response in V. virens at the initial stage

2.3 KEGG代谢途径分析

DEG被注释到91条通路中，共获得显著富集通路8条，分别涉及次生代谢产物生物合成、淀粉和蔗糖代谢、糖酵解/糖异生、碳代谢、甲烷代谢、抗生素的生物合成、丙氨酸、天冬氨酸和谷氨酸代谢、氨基糖和核苷酸糖代谢（表3、图4-C）。这些结果进一步表明，稻曲病菌菌核在形成过程中具有活跃的营养物质代谢和能量代谢。

2.4 差异表达的基因家族分析

基因家族是指在进化过程中，由一个祖先基因通过基因复制产生两个或更多的拷贝，进而发生分化形成的一类具有相似结构和功能的基因，能够编码相似的蛋白质产物。本研究中DEG共注释到180个基因家族中，占DEG总数的22.7%（图5-A），其中111个基因家族表现为显著上调，占总数的 61.67%；11个上调基因被注释到MFS转运蛋白基因家族中，富集程度最高，其他显著上调表达基因家族如糖转运蛋白（8个）、锌指转录因子（6个）、磷酸激酶（4个）等也有着较多成员受低温诱导（图5-B—E）；显著下调基因家族占比38.33%，主要包括真菌特异性转运蛋白基因家族（4个）、磷酸转移酶基因家族（3个）等（图 6）。推测显著差异数目较多的基因家族与低温诱导稻曲病菌菌核形成有关，可能参与低温胁迫的形态学发育和生理代谢响应。

2.5 与低温诱导菌核形成相关的差异表达基因

为研究稻曲病菌在低温处理后，菌核形成过程中基因表达差异，根据稻曲病菌基因组搜索获得基因编号，主要对氧化应激反应、信号转导、跨膜运输、细胞形态、生物合成等相关基因进行分析。并根据转录组数据的 FPKM 值绘制热图表示基因的相对表达量（图 7）。这些基因在不同样本间的表达量具有显著差异。

表3 差异表达基因在KEGG代谢通路中显著富集分析Table 3 KEGG enrichment significant analysis of DEGs

图5 稻曲病菌低温诱导相关基因家族的表达分析Fig. 5 Expression analysis of gene family involved in low temperature response in V. virens at the initial stage

图6 低温诱导稻曲病菌下调表达的基因家族分析Fig. 6 Expression analysis of down-regulated gene family in V. virens under low temperature

2.5.1 氧化应激反应相关基因分析 氧化应激反应能够诱导机体在遭受胁迫时产生的 ROS（O2-和H2O2等），影响丝状真菌齐整小核菌（Sclerotium rolfsii）[18]和核盘菌（Sclerotinia sclerotiorum）[19]菌核的分化过程。而NADPH氧化酶（NADPH oxidase，Nox）、超氧化物歧化酶（superoxide dismutase，SOD）、过氧化氢酶（catalase，CAT）是影响真核生物细胞ROS水平和H2O2积累的关键酶类。对低温处理后的稻曲病菌中Nox、SOD和CAT基因表达变化进行研究，发现NADPH氧化酶基因（KDB12731.1、KDB18142.1）、SOD基因（KDB13992.1、KDB18714.1）、CAT基因（KDB11441.1、KDB15626.1）在低温处理的稻曲病菌中表达量明显增加（图 7）；结合转录组数据GO途径中差异基因在氧化还原酶活性和过氧化物酶的富集，推测由于低温胁迫，稻曲病菌的氧化应激反应强烈，内部活性氧水平升高，进而诱导了菌核的分化。

图7 样品TL911和FH1016中差异显著表达基因热图Fig. 7 Heatmap of DEGs in TL911 and FH1016

2.5.2 信号转导相关基因分析 丝裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）是真核类生物信号传递通路（MAPK通路）中的关键酶，能够调控真核细胞形态和生长。转录组数据中，该类基因（KDB16903.1）在低温诱导菌核的过程中显著差异表达（图 7），表明信号传递通路被激活，能够将感知到的外界或胞内信号向下游传递，引发一系列的级联反应。

2.5.3 跨膜运输相关基因分析 主要协助转运蛋白超家族（MFS），是目前真菌主要的跨膜运输蛋白家族之一，在真菌营养吸收中起到重要作用。通过对DEG分析，发现有22个基因表达被注释到MFS蛋白家族中，结合图5-B、5-C看出MFS蛋白家族DEG在上调表达的基因家族分类中富集程度最高；其中 3个 MFS 基因 KDB15983.1、KDB18665.1、KDB15982.1上调表达量分别达到了2.46、2.14和2.11倍（图7）。另外，糖转运蛋白家族中绝大部分基因也是表现为显著上调富集（图 5）。磷脂是微生物细胞膜的重要组成部分，而磷脂酰丝氨酸合成酶（phosphatidylserine synthase，Pss）能够催化磷脂酰丝氨酸的合成，是合成细胞膜成分的关键，该基因KDB14430.1在低温处理时表达量较高（图7）。

2.5.4 细胞形态相关基因分析 Ras和 cytochrome P450（CYP450）蛋白家族能够调控菌丝的极性生长和形态发育[20-21]，稻曲病菌转录组中注释到Ras蛋白家族中的 3个基因 KDB12837.1、KDB17934.1、KDB11781.1被诱导上调表达；此外，有22个基因表达被注释到CYP450蛋白中，其中18个家族基因表达被显著诱导，以KDB17652.1表达量上调幅度最大，为7.13倍（图7）。

2.5.5 生物合成基因相关分析 真菌菌核的发育形成往往伴随黑色素的沉积，聚酮合酶是真菌黑色素合成途径中的关键酶之一，通过对 DEG进行分析，发现聚酮合酶基因（KDB13139.1）上调表达了1.65倍；此外几丁质和 1,3-β-葡聚糖作为真菌细胞壁的组成成分，其几丁质合成酶（KDB15245.1、KDB11225.1、KDB11215.1、KDB14486.1）表达量也显著增高。分析富集在淀粉和蔗糖代谢途径中的DEG，发现参与合成 1,3-β-葡聚糖的 2个关键基因（KDB12252.1、KDB16140.1）在低温处理的稻曲球中显著上调（图7）。说明代谢产物黑色素、几丁质和多糖随菌核发育而被大量合成。

2.5.6 生殖相关基因分析 菌核作为无性繁殖体是稻曲病菌分化形成的一团紧密交织的菌丝体，越冬后的菌核可进行有性生殖，萌发形成子实体并产生子囊孢子。而Velvet家族蛋白被认为是真菌特有的无性和有性发育的重要调控因子，同时也影响着次级产物的合成[22-23]。通过与模式真菌比对，找出稻曲转录组中关于Velvet family的同源基因VEA（KDB14625.1），并发现其表达水平上调了1.31倍；其中，低温处理组FPKM为36.97、对照组FPKM为28.19，明显低于整个转录组数据的平均数值130，说明VEA在稻曲病菌的基因组内属于低表达的组成型基因。另外，稻曲病菌的有性配合主要受基因MAT1-1-1和MAT1-2-1调控，发现MAT1-2-1表达量极低（FPKM value＜2），而交配型基因MAT1-1-1在低温处理组和对照组中均未检测到该基因的表达。

2.6 qRT-PCR验证基因表达

为验证转录组数据是否准确、可信，随机选取 6个DEG进行qRT-PCR检测。qRT-PCR测定结果与转录组测序输出数据分析得到的变化趋势一致（图8），表明本研究中基于转录组数据得到的DEG是可信的。

图8 随机挑选6个差异基因的实时荧光定量分析Fig. 8 qRT-PCR analysis of 6 randomly selected DEGs

3 讨论

氧化应激反应特别是ROS的形成，是诱发菌核分化和调节真菌发育的重要因素之一，高氧化状态可以驱动真菌菌核、有性结构的形成[24-26]。Nox能在真核生物中产生活性氧，是超氧化物的一个重要来源[27]。有研究表明，bcnoxa和bcnoxb参与了灰霉病菌（Botrytic cinerea）菌核的形成过程[28]。另一方面，在高氧化应激状态下，各种抗氧化酶和其他小分子（如SOD、CAT等）可能在保护真菌细胞免受ROS的有害影响方面起重要作用。SOD是细胞内ROS的第一道防线，是氧化应激过程中清除自由基最重要的酶之一，在各种生物系统中都会遇到ROS，其作用是消除O2-。CAT能够催化H2O2产生H2O和O2-，发挥抗氧化作用[29]。在本研究转录组数据中发现，CAT和SOD基因表达量升高，这可能归因于氧化还原反应的平衡被破坏和ROS含量相应地急剧增加。因此，SOD和CAT能在应激时主动清除Nox产生的活性氧，最终导致抗氧化酶浓度升高，这些结果间接证明了菌核形成与氧化应激密切相关。GRINTZALIS等研究表明，黄曲霉（Aspergillus flavus）能够在氧化应激反应中合成次级代谢产物黄曲霉素B1[30]。推测氧化应激反应不仅是菌核形成所需要的，也与菌核形成过程中次级代谢产物的合成有关联，这与本研究KEGG通路中差异基因在次级代谢产物合成途径中显著富集相吻合。

此外，MAPK信号转导通路中的丝裂霉素活化蛋白激酶基因在低温处理组中上调表达，使得ROS能够通过MAPK途径将信号传递给下游，进一步影响稻曲病菌生理过程基因的表达。例如Pss基因表达量升高，稻曲病菌通过加速合成磷脂来保护细胞膜系统并确保膜上的功能酶能够正常行使功能，结合基因家族分析，MFS转运蛋白家族和糖转运蛋白家族的DEG数目较多且表达量变化明显，说明稻曲病菌在菌核形成过程中细胞营养代谢旺盛，自身能够在膜转运蛋白的协助下实现物质的快速转运，使机体应答低温胁迫。在细胞形态发育方面，Ras是一类调控细胞极性、周期和功能分化的信号分子，而CYP450蛋白在植物木质素中间体、甾醇、萜烯、类黄酮等多种次生代谢物的合成中起关键作用[31]，也能够通过对激素的响应进而影响形态发育。Ras和CYP450的表达上调，推测其参与调控稻曲菌丝的细胞极性生长与形态发育。而锌指蛋白是一类具有手指状结构域的转录因子，在转录和翻译水平上调控基因的表达和细胞分化，如构巢曲霉（Aspergillus nidulans）中，nsdC蛋白是一种C2H2型锌指蛋白，基因敲除突变体的生长速率降低，并且完全不能形成子实体[32]。本研究锌指转录因子中6个显著上调，2个显著下调，锌指转录因子家族中的不同TF表达上调或下调，推测这些TF的种类及表达水平的变化可能与稻曲病菌对低温胁迫响应的形态发育有关。此外，代谢产物黑色素、几丁质和多糖也随菌核的发育而被大量合成。综上，稻曲病菌菌丝形成菌核无疑是对外界刺激的一种反应，通过自身形态学及生理代谢等改变对低温胁迫作出响应。

另外，在大部分丝状真菌中有一类velvet，VEA是Velvet蛋白复合物的核心蛋白，具有广泛存在性和功能保守性，能够调节菌核的形成。在寄生曲霉（Aspergillus parasiticus）和黄曲霉中VEA不仅能够调控菌核的发育、诱导黄曲霉毒素的产生，还进一步证实是机体发育和次级代谢之间分子联系的桥梁[33-34]。在烟曲霉（Aspergillus fumigatus）中VEA可以激活有性生殖，减少分生孢子的产生。但是其组成型的基因表达，使得VEA几乎不因无性或有性生殖的发生而有很大的表达变化[35]，这与本研究在 DEG分析中发现VEA较低的上调倍数相一致，故推测稻曲病菌的VEA可能以其他独特的方式调控菌核的形成。进一步挖掘稻曲病菌有性生殖基因的表达，发现在低温处理组与对照组中不仅MAT1-2-1表达量极低而且还未检测到MAT1-1-1的表达，推测是由于在低温环境下促进的是稻曲病菌菌核的快速形成，而菌核越冬后需要在适合的温度以及充足的光照下才能进行有性生殖形成子实体，所以导致在低温诱导菌核形成的过程中并未发现有性生殖基因的大量表达。由于低温诱导稻曲病菌菌核形成所涉及的分子机制尚未确定，因此通过转录组测序分析的大量基因将作为进一步研究的候选者，为今后稻曲病的防治提供数据和理论基础。

4 结论

通过对低温诱导稻曲病菌菌核形成的转录组数据进行深度挖掘，发现稻曲菌丝低温胁迫下，致使菌丝机体内部处于氧化应激状态，并通过信号转导途径放大，该过程由多个基因参与并调控基因家族成员，最终促使跨膜运输、细胞形态、生物合成等基因的上调表达，使得在形成菌核过程中蛋白表达活跃，达到合成细胞及物质的高峰期，进而促进菌核形成。