花魔芋响应软腐病菌侵染初期的转录组分析

何 斐冯明月王菁轩邵永春王 甜王 莹卢谢敏方姝颖

(安康学院 现代农业与生物科技学院,陕西安康 725000)

魔芋系天南星科魔芋属(Amorphophallus)多年生草本植物,它是目前为止自然界唯一能大量提供葡甘聚糖的重要经济作物[1]。这种天然高分子优质膳食纤维-葡甘聚糖,具有凝胶性、增稠性、持水性和成膜性等独特的特性,被广泛应用于生物医药、食品和化工等领域[2]。由于其较高的经济价值,近年来魔芋已发展成为云南、贵州、四川、陕南及湖北等山区农民脱贫致富的支柱产业[3]。花魔芋是中国的主栽品种,其产量大、葡甘聚糖含量高,但抗病性较差[4]。其中,土传和种传细菌性软腐病是花魔芋生产过程中危害最大且极难防控的病害,在植株生长期和采后贮藏期均可发病,可危害花魔芋叶片、叶柄及地下块茎,轻则造成减产20%～30%,重则导致80%以上减产甚至绝收[4],严重制约魔芋产业的可持续发展。随着科技的不断进步,利用抗病基因进行分子育种已成为解决作物抗病性的有效途径,这一途径的关键是对花魔芋新型内源抗病基因的挖掘和筛选以及进一步对其抗病防御机制的解析。

目前,对花魔芋软腐病相关的研究虽然已有一定报道,但是主要集中在病原菌鉴定与检测[3]、病害防治[5]及传统的抗病品种选育[6]等方面。Wu等[7]研究发现,芽孢杆菌产生的酰基高丝氨酸内酯酶(Aii A)通过干扰其群体感应系统的信号分子N-酰基高丝氨酸内酯(AHL),从而对花魔芋软腐病菌产生抑制作用。雷珍珍等[8]从抗软腐病花魔芋植株基因组中分离获得一条NBSLRR 类抗病基因同源序列。上述结果为花魔芋对软腐病菌侵染响应分子机制的研究提供了理论基础,但这些研究仅局限于病原菌群体感应的信号分子或花魔芋中部分抗病基因的分析,而植物对病原菌侵染的抗性反应涉及多种途径和基因的表达调控,有必要系统鉴定花魔芋中参与软腐病胁迫响应的多种基因,以揭示花魔芋抗病的分子机制及挖掘重要的抗病基因与调控功能。

植物染病后,基因转录和翻译水平通常发生明显改变,而基因表达水平的变化则早于病症的出现。李玲[9]在研究2个稻瘟病菌生理小种侵染时,发现在3个水稻品种上接种不同病菌后24 h和48 h分别有12 456个和4 851个基因发生差异表达,说明在染病早期草本植物的基因表达已经明显改变。此外,软腐病对花魔芋造成的危害严重且不可逆,增强植株自身染病初期的抗性十分必要。转录组测序技术(RNA-seq)作为一种高效快捷、准确度高、成本低的研究手段,被广泛应用于模式植物拟南芥[10]和玉米[11]、非模式植物甘蔗[12]和花生[13]等多种植物抗病机理研究。本研究通过转录组测序技术探究花魔芋感染软腐病菌初期基因表达的整体变化情况,通过对差异表达基因进行GO 功能分类和KEGG 富集性分析,筛选花魔芋响应细菌性软腐病菌侵染的相关基因、抗病相关的信号传导途径和生物活性物代谢通路,以期为深入研究魔芋抗病分子机制提供理论参考,也为魔芋抗病分子育种提供候选基因。

1 材料与方法

1.1 材料

植物材料为1 a生的安康市本地花魔芋盆栽苗。花魔芋软腐病菌菊果胶杆菌(Pectobacterium chrysanthemi)CZS-B6菌株由中国典型培养物保藏中心(China Center for Type Culture Collection,简称CCTCC,也称武汉大学保藏中心)提供。

1.2 方法

1.2.1 病原菌接种与取样 花魔芋软腐病菌的培养和接种体的制备参考徐炜等[14]的方法。参照何斐等[15]的方法进行盆栽花魔芋的种植与管理,于盆栽花魔芋生长至苗高约10 cm 左右,采用浇灌法[14]将稀释至108CFU/m L 的接种体菌悬液浇灌到健康的盆栽花魔芋植株根际,每株浇灌200 m L,以浇灌无菌水的花魔芋作为对照。根据预试验效果,分别选取浇灌处理第5天的健康和软腐病菌侵染初期的花魔芋球茎,于液氮中迅速冷冻后至-80℃保存。每6个生物学重复当作1个处理,每处理重复2次。

1.2.2 RNA 提取、cDNA 文库构建和转录组测序 将冻存样品在液氮中研磨,用RNApure Plant Kit(含DNaseⅠ)试剂盒(康为世纪生物科技有限公司,北京)提取总RNA,分别采用Nanodrop 2000 和Agilent 2100 生物分析仪(Agilent Technology,美国)检测RNA 样品的浓度和质量。经检测合格后送至北京组学生物科技有限公司构建cDNA 文库,采用实时定量PCR 方法库检合格后,使用Illumina HiSeq 2500 测序平台(Illumina Inc.,美国)进行转录组测序(RNAseq)。

1.2.3 测序数据处理 花魔芋无参考基因组,采用Trinity软件[16]对原始测序数据(Raw data)过滤后获得高质量数据(Clean data),对其进行组装拼接成转录本。去冗余后获得unigenes,使用COG、GO[17]、KEGG[18]、KOG[19]、Pfam[20]、Swiss-Prot[21]和Nr数据库进行基因功能注释。

差异表达基因筛选:用FPKM 值(每千个碱基的转录每百万映射读取的片段数)[22]计算基因表达量。差异表达基因筛选标准:花魔芋健康株与软腐病株基因表达水平差异倍数FC(Fold change)≥2,且错误发现率FDR(False discovery rate)<0.01[23]。

差异表达基因GO 注释和KEGG 富集分析:利用Blast2GO(https://www.blast2go.org/)进行差异表达基因的GO 注释[24]。通过P值(P<0.05)的Benjamini-Hochberg矫正方法确定差异表达基因中显著富集的GO 分类,使用Cytoscape(http://www.cytoscape.org/)注释通路。利用KEGG 数据库(http://www.genome.jp /kegg/),以P<0.05且Q≤0.05作为显著富集标准,进行差异表达基因的通路富集分析[25]。

2 结果与分析

2.1 测序数据质量分析

利用Illumina测序平台对健康和软腐病菌侵染初期的花魔芋球茎转录组进行测序,经过对测序原始数据质量控制,在健康和发病魔芋中分别获得19 236 352和21 116 175条clean reads(NCBI SRA 数据库登录号:PRJNA738708,各样本序列登录号 SAMN19759221 - SAMN 19759224),GC 含量分别为51.32%和51.36%,Q30值均高于93%。各样品的clean reads与拼接得到的转录本参考序列的比对效率达到75%以上,数据显示转录组测序质量良好,可进行下一步分析(表1)。

表1 花魔芋球茎转录组测序数据Table 1 Transcriptome data assembly for corm of A.konjac

2.2 Unigene功能注释

将unigene序列与COG、GO、KEGG、KOG、Pfam、Swissprot和Nr数据库比对和注释,共获得57 525条有注释信息的unigenes,占全部unigenes的51.6%。其中COG 数据库注释信息最少,有18 357条;Nr数据库注释信息最多,高达53 611条,占总注释unigenes的93.2%(表2)。

表2 Unigene注释分析Table 2 Unigene annotation based on different databases

2.3 差异表达基因分析

转录组分析(图1)表明,软腐病菌侵染初期的花魔芋球茎中共出现差异表达基因3 222个。与健株相比,病株中有2 660 个基因表达上调,562个基因表达下调,上调基因比下调表达的基因多373.3%。差异基因表达量上调、下调最大的前10个基因中(表3),多数基因编码特定功能的蛋白质或酶,如上调表达基因TRINITY_DN674_c0_g3 编码Mi AMP1 抗菌蛋白,下调表达基因TRINITY_DN3252_c0_g3 编码热激蛋白,TRINITY_DN11231_c0_g1编码胰蛋白酶与蛋白酶抑制剂、胰凝乳蛋白酶抑制剂。研究表明这些蛋白质、酶与抵御病菌侵染密切相关[26-27]。

表3 花魔芋在软腐病菌胁迫下表达量上升和下降最大的前10个基因Table 3 Top 10 up-and down-regulated genes differentially expressed in A.konjac under P.chrysanthemi stress

图1 软腐病菌侵染初期花魔芋球茎中的差异表达基因分析Fig.1 Differential expression analysis of genes in corm of A.konjac with early response to P.chrysanthemi infection

2.4 差异表达基因功能分类和富集分析

对差异表达基因进行GO 功能分类(图2),共分为3大类53个功能组。在生物学过程、细胞组分和分子功能3大类中,2个差异表达最多的功能组分别是新陈代谢过程(706个差异表达基因)和细胞进程(620个)、细胞组分(627个)和细胞(626个)、催化活性(644个)和结合(589个)。

图2 基于转录组测序数据的差异表达基因GO 功能分类Fig.2 GO functional classifications of differentially expressed genes(DEGs)based on transcriptome sequencing data

通过KEGG 对花魔芋各通路进行注释分析(图3),发现1 148个差异表达基因在细胞过程、环境信息处理、基因信息处理、代谢和有机系统5大类的信号转导及植物-病原菌互作等通路中富集。在29 条通路中,与代谢有关的通路最多(17条),占比58.6%。

图3 差异表达基因的KEGG 通路分析Fig.3 KEGG pathway enrichment analysis of differentially expressed genes

在显著富集的7条通路(表4)中,有3条与生物合成相关,2条与代谢有关,剩余的2条分别是糖酵解/糖异生和氧化磷酸化通路。富集程度最大的2 条通路分别是类黄酮生物合成[ko00941]和苯丙烷生物合成通路[ko00940]。在对花魔芋健株与软腐病株差异表达基因分类和途径分析的基础上,推测软腐病菌侵染花魔芋诱导的代谢途径中的差异表达可能与苯丙烷和类黄酮生物合成等抗性有关,对防御病菌侵染发挥重要作用[28]。

表4 差异表达基因中显著富集的通路Table 4 Significantly enriched pathways of differentially expressed genes(DEGs)

2.5 趋势表达基因通路分析

为解析花魔芋抵御软腐病菌侵染的关键通路和基因,逐条分析差异表达基因映射在KEGG 数据库中的参考通路,结果发现一些通路的差异基因存在特定的表达变化趋势(表5)。其中,硒化合物代谢和维生素B6代谢通路中,所有差异基因相比健株均呈现表达量上升的趋势。在氨基糖和核苷酸糖代谢,甘氨酸、丝氨酸和苏氨酸代谢,类黄酮生物合成等14条通路中,60%以上的差异基因相比健株表达量上调。但甾类激素生物合成(57.14%)和类胡萝卜素生物合成(66.67%)通路中的大部分差异表达基因相比健株呈现下调趋势。

表5 趋势表达基因KEGG 通路富集分析Table 5 Enrichment analysis of down-and up-regulated genes based on KEGG database

3 讨论与结论

本研究对健康和软腐病菌侵染初期的花魔芋球茎进行转录组测序分析,获得111 404条unigenes和3 222个差异表达基因,其中2 660个基因(82.6%)上调表达,562个基因(17.4%)下调表达(病株/健株)。刘梦洋[29]研究发现,大白菜抗软腐病突变体受软腐病菌侵染初期12 h 时,412 个基因(81.9%)表达量上调,91 个基因(18.1%)表达量下调。于秋佟[30]通过转录组测序发现,软腐病菌侵染24 h的拟南芥有579个基因(89.0%)表达量上调,72个基因(11.1%)表达量下调。这些研究与本结果类似,即上调表达基因数量明显多于下调表达基因的数量,说明更多的基因被诱导表达来抵御软腐病菌的入侵和发展。本试验结果为后续鉴定花魔芋抵御软腐病菌侵染初期发挥关键作用的基因缩小了筛选范围,也为魔芋抗病育种奠定基础。

在分析软腐病菌胁迫下花魔芋中差异表达变化较大的基因时(表3),鉴定到部分响应软腐病菌胁迫的差异表达基因。例如,在上调表达的基因中TRINITY_DN674_c0_g3编码MiAMP1抗菌蛋白,已有研究表明其与植物抗病有关,Mi-AMP1抗菌蛋白重组毕赤酵母能诱导梨对扩展青霉产生抗性[31]。TRINITY_DN 1022_c0_g3编码类枯草杆菌蛋白酶家族,有研究显示类枯草杆菌蛋白酶在病原物侵染后被特异性诱导并参与宿主植物细胞程序性死亡,同时在病原物识别和信号转导级联反应中发挥重要作用[32]。值得关注的是,笔者发现也存在与抗病相关的下调表达基因,如TRINITY_DN3252_c0_g3编码热激蛋白,已有报道显示众多热激蛋白参与调控植物抗逆和防卫反应[33]。TRINITY_DN11231_c0_g1 编码胰蛋白酶和蛋白酶抑制剂,多种植物的蛋白酶抑制剂具有抗病虫功能[34]。据此,推测花魔芋可能根据病菌胁迫的特点,通过上调表达具有针对性的抗性基因,同时下调表达广谱抗性基因来响应软腐病菌的胁迫。

本研究在分析KEGG 代谢通路时,发现花魔芋中16条通路的差异表达基因在软腐病菌胁迫下全部或大部分上调表达(表5)。经逐条分析,可按照功能将其分成两类:第1类与宿主植物抗逆有关;第2类与细胞合成原料相关。在第1类通路中,硒化合物和维生素B6代谢通路中差异表达基因全部上调。已有研究证明,硒是植物生长发育不可缺少的微量元素,其通过维持植物体内较高的抗氧化酶活性来诱导宿主产生抗病性[35];维生素B6作为多种代谢酶的重要辅助因子,参与植物氨基酸、乙烯合成等多种基本代谢活动,亦能通过调整细胞内抗氧化能力从而调控番茄植株对灰霉病的抗性[36]。黄酮类化合物是植物自身抗病、抗逆的主要成分[37],在本研究中,涉及黄酮类化合物合成的类黄酮生物合成通路中有91.0%差异基因上调表达。N-聚糖生物合成通路合成的寡糖是一类植物抗逆诱导剂[38],其中木寡糖能诱导拟南芥对半活体病原细菌丁香假单胞菌番茄致病变种、活体病原烟草花叶病毒和死体病原真菌核盘菌产生抗性[39]。在本研究中N-聚糖生物合成通路83.3%差异基因上调表达。而链霉素生物合成通路合成的链霉素主要用于防治由革兰氏阴性细菌引起的烟草野火病和青枯病、大白菜软腐病、柑橘溃疡病、水稻白叶枯病和条斑病等植物病害[40],在本研究中链霉素生物合成通路60%差异基因上调表达,推测该通路合成的链霉素可能在抵御软腐病菌侵入过程中发挥重要作用。第1大类参与宿主植物抗逆相关的信号转导和物质运输的主要通路还包括氧化磷酸化、MAPK 信号通路、氨基糖和核苷酸糖代谢通路、植物-病原菌互作和植物激素信号转导[41],这些通路75%以上的差异基因相比于健株表达量上调。

在第2类通路中,半胱氨酸和蛋氨酸代谢,甘氨酸、丝氨酸和苏氨酸代谢,谷胱甘肽代谢,苯丙氨酸代谢,苯丙氨酸、酪氨酸和色氨酸的生物合成均与氨基酸的代谢或合成相关,这些通路均有超过66.7%的差异基因受诱导上调表达,这可能引起相应的氨基酸合成加快,从而为合成抗病相关的蛋白和酶提供充足原料[42]。然而甾类激素和类胡萝卜素生物合成通路的大部分差异基因表达下调(表5)。广泛存在于植物体内的类胡萝卜素和甾类激素主要参与调控植物生长发育及生理功能[43],且类胡萝卜素含量也与植物抗病性相关[44]。类胡萝卜素和甾类激素生物合成通路中基因表达的下调可能引起细胞内甾类激素和类胡萝卜素供应紧缺,从而影响细胞包括抗软腐病相关基因表达等在内的正常生理功能,这可能是染病花魔芋无法有效抵御软腐病菌入侵的原因之一。

本研究通过RNA-seq测序得到3222个在花魔芋感染细菌性软腐病初期出现差异表达的基因,筛选出与抗病相关且表达量差异大的基因。进一步进行GO 注释和KEGG 通路分析发现2条通路的富集程度最大且合成物质具有抗病性-苯丙烷生物合成和类黄酮生物合成通路。通过这2条通路合成的物质可能在花魔芋抵御软腐病菌侵入过程中发挥关键作用。同时,研究发现16条通路中60%以上的差异基因上调表达,而2条通路中大部分差异基因下调表达。其中,硒化合物代谢、维生素B6代谢、类黄酮生物合成、N-聚糖生物合成及链霉素生物合成这5条通路生成的相关物质与花魔芋抗软腐病菌胁迫密切相关。这些抗病相关通路和差异表达基因的分析结果为花魔芋抗病机制的研究奠定了基础,也为抗病相关基因的功能研究提供一定参考。