油茶炭疽病菌果生刺盘孢效应子的筛选*

陈星州 周国英 陈行钢 江玲玉 包安华 刘君昂

(中南林业科技大学 南方人工林病虫害防治国家林业局重点实验室 森林有害生物防控湖南省重点实验室经济林培育与保护教育部重点实验室 长沙 410004)

油茶(Camelliaoleifera)作为农村兴林富民的主导产业和脱贫攻坚的民生产业，在我国种植面积不断扩大。炭疽病是油茶栽培区主要病害之一，引起油茶落果、落蕾，甚至整株衰亡，在湖南省油茶种植区引起的落果率平均为20%，有时可高达40%～50%(喻锦秀等， 2014)。前期调查研究发现，中国油茶炭疽病的病原菌种类复杂多样，果生刺盘孢(Colletotrichumfructicola)是目前各产区油茶炭疽病的优势病原菌(李河， 2018)。

自然条件下，植物体会受到多种病原体的攻击，因此植物形成了复杂的免疫系统，植物和病原体之间的相互作用现在被描述为一种共同进化的军备竞赛，被称之为之字形模型(Jonesetal.,2006)。当病原体攻击植物细胞时首先会促发植物一般防御机制，这个过程被指定为PAMP(pathogen-associated molecular pattern)触发免疫模式或者触发免疫(PAMP-triggered immunity,PTI)(Kushalappaetal.,2016)。另外，病原体也会通过将效应子分泌到宿主细胞中来干扰PTI，植物通过直接或间接识别病原体效应物的抗性(R)蛋白激活第二层免疫系统，这被称为效应物触发免疫(effector-triggered immunity,ETI)。病原体通常会进化出新的效应子，并转移到宿主细胞中干扰植物的免疫系统，从而促进病原体定殖(Dangletal., 2013； Yanetal., 2018)。例如果生刺盘孢在早期侵染苹果(Maluspumila)树时分泌的效应子CfEC96，缺失该基因的果生刺盘孢突变体的营养生长、附着胞的形成均不受影响，但是突变体对苹果叶和果实的毒力减弱(Shangetal., 2020)。真菌效应子都具有相似的分子特征:无糖基磷脂酰肌醇(Glycosylphosphatidyl inositol, GPI)的锚定位点，无跨膜结构域，没有将蛋白输送到线粒体或其他胞内细胞器的预测定位信号，含有N-端信号肽，富含半胱氨酸，氨基酸残基<300个，序列高度特异(Ramachandranetal., 2017)。对于真菌效应子的研究目前主要集中于稻瘟菌(Magnaportheoryzea)、番茄叶霉菌(Cladosporiumfulvum)、尖孢镰刀菌(Fusariumoxysporum)等，但对于炭疽菌效应子的研究报道较少。随着计算机技术的高速发展以及测序技术的不断更新，利用生物信息学手段对效应子进行预测与分析逐渐成为筛选真菌候选效应子的重要方法，并且有证据证明转录组学的数据可以预测毒力因子(Jonesetal., 2018)。

本研究综合利用BUSCA、Target P、big-PI Predictor、GO、KEGG功能富集和PHI等生物信息分析软件和数据库，对果生刺盘孢分泌蛋白和候选效应子进行预测和分析，并通过烟草(Nicotianabenthamiana)瞬时表达系统验证候选效应子的功能，拟为果生刺盘孢致病相关基因研究提供理论依据，同时为油茶种植的病害管理提供参考。

1 材料与方法

1.1 试验材料

野生型本氏烟草(Nicotianabenthamia)种子、果生刺盘孢菌、质粒pEGAD载体由中南林业科技大学森林保护重点实验室保存。大肠杆菌(Escherichiacoli)DH5α、农杆菌(Agrobacteriumtumefaciens)GV3101由南京诺唯赞生物科技股份有限公司提供。

1.2 果生刺盘孢转录组数据来源

转录组测序由广州基迪奥生物科技有限公司完成，测序样品包括果生刺盘孢分生孢子组织以及由分生孢子侵染油茶叶片96 h后的病斑组织，每组样品设置3个生物学重复。参考基因组为ColletotrichumfructicolaNara gc5基因组(BioProject Accession： PRJNA225509)。差异表达基因的筛选阈值为相对表达量|log2(fc)|≥1且FDR≤0.05，最终得到7 850个侵染前后差异表达的基于蛋白氨基酸序列用于后续分析。

1.3 果生刺盘孢经典分泌蛋白的生物信息学分析

将同时含有N端信号肽、亚细胞定位为胞外分泌类型、无跨膜结构域、无GPI锚定位点4个特征的蛋白质定义为果生刺盘孢经典分泌蛋白。

BUSCA(http∥busca.biocomp.unibo.it)(Castrenseetal., 2018)可同时对果生刺盘孢转录组中编码的全部蛋白序列的信号肽、跨膜结构域、亚细胞定位、GPI锚定位点进行分析。利用TargetP 1.1 Server(http:∥www.cbs.dtu.dk/services/TargetP/)对上述定位到胞外的蛋白进行验证分析(Emanuelssonetal., 2007)。利用LipoP 1.0 Server(www.cbs.dtu.dk/services/LipoP/)对经典分泌蛋白信号肽切割位点的氨基酸组成进行分析(Junckeretal., 2003)。

1.4 果生刺盘孢经典分泌蛋白的功能注释

利用GO(http:∥www.geneontology.org)、KEGG(http:∥www.genome.jp/kegg/)(Kanehisaetal., 2019)、PHI(http:∥www.phi-base.org/)(Martinetal., 2019)多个数据库对经典分泌蛋白序列进行BLASTP比对，对其功能进行注释。

1.5 果生刺盘孢候选效应子的预测分析

对果生刺盘孢经典分泌蛋白进一步分析筛选，其筛选标准为氨基酸长度<300且带有4个或多个半胱氨酸残基(Jonesetal., 2018)，并利用Effector P(http:∥effectorp.csiro.au/)在线分析软件对上述蛋白质序列进行分析(Pengetal.,1999)。

1.6 候选效应子qRT-PCR定量分析

随机选取9个候选效应子，利用qRT-PCR技术验证候选效应子的基因表达情况(表1)。按照RNAprep Pure多糖多酚植物总RNA提取试剂盒(天根生化科技(北京)有限公司)说明，分别提取分生孢子和叶片中的菌丝RNA。按照FastKing一步法除基因组cDNA试剂盒(天根生化科技(北京)有限公司)操作步骤，进行反转录； 以反转录产物cDNA 作为模板，按照SuperRealPreMix Plus(SYBR Green)(天根生化科技(北京)有限公司)试剂盒说明书进行qRT-PCR 分析。将分生孢子中目的基因表达量定为1，通过扩增获得Ct值来计算目的基因在侵染时期的相对表达量，Actin为内参基因(李司政等， 2020)。每个基因设置4个重复。

表1 9个候选效应子基因qRT-PCR的特异性引物序列Tab.1 9 candidate effectors gene-specific primer sequences designed for qRT-PCR

1.7 候选效应子功能验证

以cDNA作为模板，从上述9个候选效应子中随机选取4个，扩增其全长基因片段。构建pEGAD-CfEP(Colletotrichumfructicolaeffector protein)表达载体(表2)。将表达载体转化进大肠杆菌感受态DH5α中，选取测序正确的菌落提取质粒，将质粒转化进农杆菌感受态内。

携带重组载体的农杆菌加入LB液体培养基中，摇床28 ℃、200 r·min培养24～48 h。6 000 r·min-1离心2 min收集菌液，弃上清液，加入适量烟草注射液[MES(10 mmol·L-1，pH5.6)、MgCl2(10 mmol·L-1)、AS(10 mmol·L-1)]悬浮农杆菌，稀释菌液至OD600值为0.5～1.0； 将悬浮液置于28 ℃，避光孵育2～3 h。使用1 mL的无菌注射器(去除针头)，吸取适量悬浮菌液，利用渗入法将其注射到条件至适宜的本氏烟草叶片的背面。注射完成后，将烟草置于25 ℃继续培养，2～3天后观察烟草叶片表型，观察叶片坏死情况。

剪下烟草叶片，利用3，3-二氨基联苯胺(3，3-diaminobenzidine,DAB)对烟草叶片进行组织染色，具体参照何燕华等(2014)的方法。

表2 4个候选效应子基因全长克隆引物序列Tab.2 4 candidate effector gene-primer sequences designed for full-length cloning

2 结果与分析

2.1 果生刺盘孢分泌蛋白的预测

2.1.1 经典分泌蛋白的预测分析 本研究利用BUSCA软件同时对果生刺盘孢转录组中编码的7 850条蛋白质氨基酸序列的信号肽、跨膜结构域、亚细胞定位、GPI锚定位点进行预测分析，结果表明，分析的7 850条序列中有436个氨基酸序列含有N端信号肽，所占比例为5.55%； 其中含有分泌型信号肽的蛋白质中有390个蛋白不具有跨膜结构域； 在没有跨膜结构螺旋的蛋白质中有354个蛋白序列不具有GPI锚定位点； 具有以上3个特征的354个蛋白质亚细胞定位均为胞外分泌类型。进一步用TargetP 1.1 Server对分泌到胞外的蛋白质进行验证分析可知，含有信号肽的蛋白质345个，所占比例为97.46%。综合上述分析，在果生刺盘孢侵染阶段的转录组所编码的7 850条蛋白质氨基酸序列中，有345个符合经典分泌蛋白特征，所占比例为4.39%。

2.1.2 经典分泌蛋白的特征分析 1) 氨基酸长度 经典分泌蛋白氨基酸长度在各区间段分布较为均匀，而符合效应子筛选条件(<300 aa)的蛋白约占蛋白总数的36%，数量相对较少(图1A)。

2) 信号肽长度 在345个经典分泌蛋白中，信号肽长度(氨基酸残基数)最大为33个氨基酸，最小的为15个氨基酸； 其中，有64个蛋白含有长度为18个氨基酸的信号肽，为数量最多(图1B)。

图1 经典分泌蛋白的特征分析Fig.1 Characteristic analysis of classic secreted proteinsA： 氨基酸长度分析； B： 信号肽长度分析。A： Amino acid length analysis； B： Signal peptide length analysis.

3) 信号肽切割位点氨基酸出现频率及分布 根据信号肽酶识别位点的不同，可将信号肽分为信号肽酶Ⅰ型(SPⅠ)和信号肽酶Ⅱ型(SPⅡ)。SPⅠ分泌型信号肽的典型结构特征是C端切割位点(-3～-1位置)为A-X-A； SPⅡ为脂蛋白信号肽，典型结构是C端切割位点(-3～-1位置)为L-(A/S)-(A/G)。305个经典分泌蛋白含有SPⅠ型信号肽识别位点，40个含有SPⅡ型信号肽识别位点。

信号肽中氨基酸残基的分布统计结果表明，20种氨基酸的出现频率从高到低依次为ASPVLTGQIDNREFKCHYMW。其中，非极性氨基酸G、A、V、L、P的出现频率相对较频繁，A所占比例最大，为22.88%。此外，分泌蛋白在-1位的氨基酸种类较为保守，以非极性疏水氨基酸为主(A最多，所占比例高达66.09%)； 其次为-3位(P最多，所占比例为44.93%)(图2)。

图2 果生刺盘孢经典分泌蛋白中氨基酸的比例Fig.2 The ratio of amino acids in C. fructicola classic secreted proteins

2.2 果生刺盘孢经典分泌蛋白的功能预测

2.2.1 GO功能注释 利用GO数据库对345个经典分泌蛋白进行功能注释分析，共得到465条GO术语信息； 其中分子功能(molecular function)>生物途径(biological process)>细胞组分(cellular component)。在分子功能中最多的是催化活性(catalytic activity)，占22.18%(图3A)； 在生物途径中，最多的是初级代谢过程(primary metabolic process)，占25.0%(图3B)； 在细胞组分中，最多的是细胞内膜结合细胞器(intracellular membrane-bounded organelle)，占37.93%(图3C)。各组分中占比较多的蛋白质功能均与真菌的致病力以及生长发育有着直接或间接的联系。

2.2.2 KEGG功能富集 利用KEGG对经典分泌蛋白345条氨基酸序列进行分析，最终得到164条Pathway信息。其中最多的是代谢途径(metabolic pathways)占26.2%，其次是次生代谢产物(biosynthesis of secondary metabolites)及淀粉和蔗糖代谢(starch and sucrose metabolism)分别占10.98%和6.10%(图4)。在经典分泌蛋白中，碳水化合物代谢、抗氧化通路及水解酶基因等明显富集。有研究表明侵染过程中果生刺盘孢行使包括“增长”“细胞”“对刺激的响应”“代谢过程”和“信号传递”功能的蛋白明显增多(Keetal.,2014)，这与本研究结果相似。

2.2.3 PHI功能注释 通过利用PHI(pathogen-host interaction data-base)致病菌数据库，对已经获得的345条经典分泌蛋白进行功能注释，试找出经典分泌蛋白在致病菌中多行使何种功能。最终检索得到80条结果，其中31个比对序列通过反向遗传技术验证了其与致病的相关性，包括缺失后致死相关的基因1个，缺失后丧失致病力的基因4个，缺失后降低毒力的基因25个，还包含了1个已知效应子(表3)，该效应子发现于稻瘟病菌中，但在本研究检索中未发现已知的炭疽病效应子的信息，大多数已知致病基因同时在多个不同宿主致病菌中被发现，其中有一些蛋白在宿主为其他哺乳动物、鱼类和昆虫的致病菌中也有发现，但数量较少(3.75%)，例如鼠类的念珠菌(Candidaalbicans)病。

2.3 果生刺盘孢效应子预测分析

真菌效应子由于缺乏保守的氨基酸序列特征，预测一般采用相对宽泛的标准(Sperschneideretal., 2015)。大多数已知的真菌效应子应具有信号肽、分子量小以及富含半胱氨酸残基等特点。本文在已预测了果生刺盘孢经典分泌蛋白的基础上，以氨基酸长度<300且半胱氨酸残基的数量≥4为筛选条件，进一步对效应子进行筛选。123个分泌蛋白氨基酸长度符合筛选条件； 进一步分析发现，有65个分泌蛋白同时满足半胱氨酸残基≥4的条件； 即果生刺盘孢在侵染阶段有65个候选效应子差异表达，占转录本蛋白总数的0.8%，占经典分泌蛋白总数的18.84%。对效应子氨基酸长度分析发现，有64个效应子氨基酸长度小于200 aa，59个长度在200～300 aa之间。对效应子半胱氨酸残基数进行统计，发现其主要集中在4～8之间，其中以4个和6个较多，分别占27.7%和24.61%。对候选效应子进行功能分析，有36个候选效应子在ColletotrichumfructicolaNara gc5参考基因组中具有注释，其中包括几丁质结合蛋白、短链脱氢酶、酪氨酸酶以及含有CFEM结构域的蛋白(含有该结构域的蛋白质仅存在于真菌中)等。除此之外，另有29个为假定蛋白(hypothetical protein)。

此外，前人报道Effector P作为一类真菌效应子的在线预测软件，综合考虑了所预测蛋白的序列长度、分子量、蛋白电荷、半胱氨酸丝氨酸和色氨酸含量外，同时将真菌蛋白质在植物中的表达特征作为效应子预测的优选项，从而大幅度提高了效应子预测的准确性。利用Effector P在线预测软件对65个候选效应子分析发现，有31个符合Effector P效应子标准，而其中除了14个获得注释外，其他17个均为假定蛋白。

2.4 果生刺盘孢候选效应子的表达量分析

采用随机抽样的方法选取9个候选效应子基因进行qRT-PCR分析(图5)。9个基因的表达趋势与在RNA-seq数据中发现的趋势相似。侵染时期9个候选蛋白基因，与转录组数据结果相符，说明转录组的数据是可靠的。同时也符合效应子一般在侵染阶段会大量表达的情况。

图5 果生刺盘孢中9个候选效应子的定量分析Fig.5 Quantitative analysis of nine candidate effectors from C. fructicola transcriptomeRNA-Seq和qRT-PCR比较验证结果。误差条表示平均值的误差值。Comparison of RNA-Seq and qRT-PCR validation results.Error bars represent the standard error of the mean.

图6 果生刺盘孢中4个候选效应子的功能验证Fig.6 Functional verification of 4 candidate effectors in C. fructicola

2.5 果生刺盘孢候选效应子的功能验证

利用烟草瞬时表达系统验证选取的4个候选效应子均能诱使烟草叶片发黄或皱缩(图6)，且与阳性对照Bax蛋白造成的结果一致，而注射空载体的区域未发生相应变化。DAB染色实验注射候选效应子区域和Bax蛋白区域均有深褐色沉淀，注射空载体区域未发生变化，这一结果也表明注射区域烟草细胞出现损伤，有大量H2O2沉积。以上结果均表明，果生刺盘孢候选效应子基因成功导入烟草中，并激活烟草细胞的防御反应，诱使烟草细胞凋亡，与其他研究人员验证的大多数效应子功能相符。

3 讨论

3.1 果生刺盘孢经典分泌蛋白预测与筛选

分泌蛋白作为一类重要的致病因子，与植物病原真菌的致病性密切关联，随着越来越多病原真菌基因组测序的完成，应用生物信息手段对其分泌蛋白进行预测分析，有利于全面了解其致病机制。而随着生物信息学分析软件的不断更新，蛋白质数据库中数据的日益丰富，对于蛋白基本特征的预测手段也越来越多，越来越准确。通过分析草莓(Fragariaananassa)与果生刺盘孢的转录组数据(包括24、48、96 h 3个时间段)，并利用Signal P、TMHMM、TargetP和big-PI predictor预测工具分别对蛋白质信号肽、跨膜结构域、GPI锚定位和亚细胞定位分析最后得到52个候选效应子(Emanuelssonetal., 2000)，在分析过程中还含有多个来自其他植物病原体的已知效应子，如基因CGGC5_2199(Cmu1的同源物)，可抵消水杨酸(SA)宿主的依赖性免疫，降低SA前体分支糖酸盐的水平(Zhangetal., 2018)，基因CGGC5_10914(Ecp6的同源物)可以螯合几丁质寡糖以防止引起植物免疫力(Djameietal., 2011)，其结果与本文的分析结果相似。综上所述，本文利用BUSCA分析软件可以同时分析多种结构，在一定程度上减轻了预测工作的难度，同时预测结果可靠性较高。

3.2 果生刺盘孢经典分泌蛋白功能注释及分析

真菌在侵染寄主植物时会分泌大量的毒力因子，如效应子、植物细胞降解酶和次生代谢物生成酶等。多数情况下，并不是单一的毒力因子起作用而是多种酶类共同作用，致使真菌成功侵染寄主植物，此外除了毒力因子之外，真菌会产生许多帮助其生长增殖的酶类，以致于更好地完成侵染。为研究候选效应子在侵染阶段中的作用，本文利用GO、KEGG、PHI 3个数据库对果生刺盘孢侵染阶段转录组中345个分泌蛋白进行功能注释分析分别获得465条GO术语信息、164条Pathway信息和80条致病相关基因，与前人分析获得的转录组数据中富含小分泌蛋白、碳水化合物活性酶和次生代谢产物合成酶的结果相符(Jongeetal., 2010)； 就KEGG功能富集结果表明，大多数蛋白与真菌代谢途径有关，其次真菌的次生代谢产物的合成、内质网中的蛋白质的合成、淀粉和蔗糖代谢和酪氨酸代谢功能有明显的富集。因此，推测大多数蛋白的代谢途径都与真菌致病能力和维持真菌自身的生长发育相关，如酪氨酸酶就与部分真菌孢子形成以及稳定性、毒力和黑色素形成有密切的关系(Liangetal., 2018)； Halaouli等(2006)对胶孢炭疽菌转录组信息进行GO和KEGG功能富集分析时发现，碳水化合物代谢、抗氧化通路及水解酶基因等明显富集，这与本文结果相似。

3.3 果生刺盘孢候选效应子的定量分析及功能验证

在345个经典分泌蛋白中筛选得到17个未被注释的候选效应子，随机选取9个候选效应子在侵染阶段都有明显的表达，与转录组结果相似都上调表达，这与Zhang等(2016)研究结果相似。选取的4个候选效应子在烟草瞬时表达试验中也表现出具有诱使细胞坏死的功能，这进一步证实了本文所获得的候选效应子的可靠性。基于现有的研究基础，今后可深入研究候选效应子的功能和作用机制，包括利用亚细胞定位技术探究其作用位点、候选效应子的信号肽功能、筛选与候选效应子相匹配的靶标蛋白等。

4 结论

果生刺盘孢侵染油茶叶片的转录组数据中分析筛选得到17个候选效应子，qRT-PCR验证了其中9个候选效应子的表达情况同时验证了转录组数据的可靠性，通过烟草瞬时表达试验，鉴定出其中4个候选效应子具有诱导烟草叶片细胞坏死的功能。利用生物信息的手段设置合理的筛选条件对病原真菌侵染植物阶段的转录组数据进行预测分析，能有效筛选得到病原真菌效应子，这为植物与病原真菌的互作提供研究思路。