长足大竹象取食胁迫下慈竹WRKY转录因子的鉴定及生物信息学分析

唐 昊，郑 莉,张 羽，李沅秋，罗朝兵

(1.乐山师范学院生命科学学院，四川乐山 614000；2.竹类病虫防控与资源开发四川省重点实验室，四川乐山 614000)

0 引言

WRKY转录因子是植物中的调节蛋白中的一个大家庭，WRKY转录因子约有60个高度保守的氨基酸，其特征为N端WRKYGQK基序和C端CCHC- (C2HC)或CCHH- (C2H2)型锌指基序，WRKY蛋白通过与基因启动子中的W-BOX motif (TTGACC/T)结合来调控基因表达[1-4].同时根据WRKY结构域的数量和锌指基序的特征，可以划分为三个主要的家族[1-3]，其中家族Ⅰ具有两个典型的WRKY结构域和一个C2H2-型锌指基序，而其他大多数的WRKY蛋白只含有一个WRKY结构域，归属于家族II或家族Ⅲ，家族II的锌指基序与家族Ⅰ相同，家族Ⅲ的锌指基序为C2HC.此外，家族II中的WRKY蛋白可以进一步分为IIa 、 IIb、IIc、IId 和IIe五个亚族，而家族III中WRKY蛋白分为IIIa和IIIb亚族[3,5].大量的WRKY转录因子在多种植物的基因组中被发现并鉴定，例如拟南芥[1]，水稻[6]，玉米[7]和小麦[8]等.同时遗传和分子生物学研究表明，WRKY转录因子在植物中一方面参与抵御诸如细菌[9-11]、真菌[12-13]、病毒病原体[14]及病毒[15]各种生物胁迫的调控中起着关键作用；另一方面也在响应如干旱胁迫[8,16-18]、水胁迫[19]、盐胁迫[20]和寒冷胁迫[21-22]等非生物胁迫时的植物体内的调控.此外，WRKY转录因子还能参与植物自身的生长发育调控，如叶片的衰老调控[23-24]、开花调控[25-27]、细胞内信号转导[28]、植物种子的发育及毛状体起始和衰老[3,29].

慈竹(Bambusaemeiensis)属于丛生竹类，广泛分布于我国南方地区，被广泛应用于造纸及园林绿化等领域[30].其中在造纸领域，按照竹种对造纸质量的影响进行分级排序，竹种可以被分为四级，而慈竹属于最重要的第一级[31].长足大竹象(Cyrtotrachelusbuquet)是主要的竹林害虫之一，属于鞘翅目象甲科.成虫期的长足大竹象可将喙完全的伸入幼嫩竹笋内部, 取食竹腔内层的幼嫩组织[32]，造成新竹生长不良或形成脓包竹，引起竹腔变黑、腐烂，造成大量退笋，这将极大的影响了慈竹的产业化发展[33-34].

根据目前的研究现状，对于慈竹而言，其目前尚未有基因组草图公布，且暂未有慈竹在长足大竹象取食胁迫下的WRKY转录因子的研究.因此本文结合长足大竹象取食胁迫状态下慈竹转录组的测序及生物信息学分析，对慈竹中WRKY转录因子进行鉴定及分析，以期发现WRKY转录因子在慈竹抗虫害表征中所发挥的作用，同时为研究慈竹WRKY转录因子结构和功能提供信息参考.

1 材料与方法

1.1 植物和长足大竹象材料的准备和处理

新鲜慈竹竹笋样本采集于四川省乐山市沐川县.长足大竹象(C.buqueti)于2017年8月收集于四川省沐川县，成虫在出土后三日喂养于光照培养箱中(保持26℃恒温和70%的湿度，并维持16 h的光照及8 h的黑暗处理).在试验处理前，使长足大竹象保持饥饿状态24 h.

取100只饥饿处理后长足大竹象成虫置于新鲜竹笋上，分为长足大竹象未取食的竹笋部位(Group1)及同一竹笋上长足大竹象取食的竹笋部位(Group3)及同一丛竹笋中另一竹笋上长足大竹象完全未取食的竹笋部位(Group2)三个组，Group1和Group2取样的部位尽量保持在不同竹笋上近似的部位.实验的详细处理方法参照李沅秋等的实验方法[35].

1.2 慈竹转录组中WRKY转录因子序列的获得

提取慈竹竹笋RNA，经质量检测后用Illumina HiSeqTM2000平台测序.将测序得到的原始数据过滤掉低质量和冗余的Reads读段后，采用 Trinity 软件进行数据的从头组装及拼接，最终得到慈竹的转录本 Transcripts并从中选取最主要的转录本作为 Unigene序列，进一步并对 Unigene 数据去冗余处理并得到非冗余 Unigene 库.对Unigene序列进行基因功能注释[36,37].根据注释结果搜索“WRKY transcription factor”，查找注释相关的Unigene名称，手动调取 Transcript.fa中上述Unigene序列，去除掉序列小于200bp的Unigene基因序列信息.

1.3 慈竹WRKY转录因子家族成员的鉴定及其蛋白质的性质分析

利用ORFFinder( https://www.ncbi.nlm.nih.gov/orffinder/)功能查找每个基因的ORF，只保留氨基酸残基数大于100的相应核苷酸序列作进一步分析.利用ProtParam程序(http://web.expasy.org/protparam/)对WRKY蛋白进行相关性质参数的预测[38].

1.4 系统进化树的构建及保守motif 的分析

在Plant Transcription Factor Database(http://planttfdb.cbi.pku.edu.cn/.)中下载拟南芥(Arabidopsisthaliana)WRKY蛋白序列，使用DNAMAN6.0及MEGA 7.0软件对慈竹及山羊草的WRKY蛋白进行多序列比对，然后在MEGA 7.0软件中的系统进化树构建方法中选用邻接法并设置1000次重复bootstrap检验、Poisson模型构建WRKY蛋白的系统发生树[39-40].利用在线网站(http://meme-suite.org/index.html)并使用默认参数设置对WRKY转录因子蛋白进行motif分析[41].

1.5 慈竹WRKY表达模式

根据慈竹转录组数据中鉴定出的WRKY转录因子相对表达量的丰度值，将值导入 MeV 4.9.0分析软件中，绘制热图并进行聚类分析[42].

2 结果与分析

2.1 WRKY转录因子蛋白的性质分析

根据ORFFinder查找结果，从慈竹转录组中筛选出29个大于100个氨基酸的ORF序列.这些WRKY蛋白的ORF的氨基酸含量差别较大，最长ORF的氨基酸数从153(comp46865_seq1)到688(comp7071_seq3)；相对分子量位于16234.32 Da(comp46865_seq1)到74034.74 Da(comp7071_seq3)之间；等电点介于4.89(comp55566_seq0)到10.15(comp46865_seq1)之间；此外，蛋白质的亲水性总平均值结果表明慈竹WRKY 蛋白大多为亲水性蛋白质(表1).

表1 慈竹WRKY 蛋白的相关性质分析Tab.1 Analysis of the related properties of the WRKY protein in Bambusa emeiensis

2.2 慈竹WRKY转录因子蛋白的结构域分析

使用DNAMAN6.0软件对初步鉴定出的29条慈竹WRKY转录因子蛋白进行多序列比对，结果如下：在29条WRKY蛋白质序列中，有4条序列C端缺失(comp9372_seq1、comp17764_seq1、comp34208_seq1及comp51858_seq0)，主要原因或为转录组数据拼接不完整所致；其余25条具有完整WRKY转录因子蛋白质结构域的序列中，有8条序列C端为CCHC型锌指结构，17条为C端CCHH型锌指结构；在N端， comp34208_seq1的WRKYGQK基序变为了WPKYGSQK，存在着2个氨基酸残基的变异.此外，comp5390_seq0和comp51858_seq0转录因子的WRKYGQK基序中，Q(谷氨酰胺)被K(赖氨酸)所替代(图1).

图1 慈竹WRKY转录因子序列比对分析Fig.1 Sequence alignment of Bambusa emeiensis WRKY transcription Factor

2.3 慈竹及拟南芥WRKY转录因子的系统进化树构建

图2 慈竹与拟南芥WRKY转录因子的系统进化树Fig.2 Phylogenetic tree of WRKY Transcription Factor in Bambusa emeiensis and Arabidopsis thaliana

为进一步对慈竹WRKY转录因子的进化关系进行分析，从Plant Transcription Factor Database下载了41条拟南芥(Arabidopsisthaliana)WRKY转录因子的蛋白质序列，利用MEGA 7.0软件中CLUSTAL参数对这70条WRKY蛋白序列进行比对，绘制出慈竹WRKY 转录因子的进化树(图2).

根据系统进化树的分析结果，可以将70条序列分为3个家族.其中家族Ⅰ(Group 1)含有6个慈竹WRKY转录因子序列，家族II(Group 2)含有包括5个亚族，共有15个WRKY转录因子序列.家族III(Group 3)包含8个序列慈竹WRKY转录因子.

2.4 慈竹WRKY家族motif分析

利用MEME分析慈竹WRKY的motif(图3)，结果发现所有慈竹WRKY转录因子均含有motif1，其序列包含WRKYGQK基序，这与多序列比对结果相吻合(图1).

此外、comp34208_seq1、comp9372_seq1、comp17764_seq1缺失motif3和motif5，导致它们的C端CCHC- (C2HC)或CCHH- (C2H2)型锌指基序的缺失，可能原因是在进行序列拼接时组装不完善所致.

2.5 慈竹WRKY基因表达分析

目前尚未有慈竹WRKY转录因子在虫害(长足大竹象)取食胁迫下基因表达量变化的研究，根据本试验的结果(图4)：在不同竹笋的未被取食部位(Group1与Group2)，WRKY转录因子在慈竹中表达量基本接近，而在同一根竹笋的取食与未被取食部位(Group1与Group3)，comp5294_seq2、comp7144_seq0、comp4915_seq1、comp1871_seq0、comp5390_seq0、comp26614_seq1、comp30333_seq1等WRKY转录因子表达量存在着较大的变化，表明这几个WRKY转录因子可能对慈竹抵御虫害有重要作用.

图3 慈竹WRKY转录因子motif分析Fig.3 Analysis of Bambusa emeiensisWRKY transcription factor motif图4 长足大竹象取食胁迫下慈竹WRKY转录因子的表达模式分析Fig.4 Analysis of expression patterns of Bambusa emeiensis WRKY transcription factors under under the Cyrtotrachelus buquet feeding-stress

3 讨论与结论

WRKY因子通常被认为是存在的植物特异性，但在驯化和多倍体化的选择过程中，野生和栽培植物中都保留了重复的WRKY转录因子[43].此外，WRKY转录因子被发现与倍间杂交种的生存能力有关[44].系统发育序列分析和比较转录组学表明，它们在单子叶植物和双子叶植物之间保留了各自的功能[45].WRKY转录因子家族在多种模式植物的基因组中被发现并鉴定，如拟南芥[1]，水稻[6]和玉米[7]等.在慈竹中，刘红梅等鉴定了两个 WRKY转录因子，它们主要参与慈竹对逆境的响应及干旱和ABA(脱落酸)的胁迫响应[46],而对于慈竹中WRKY转录因子在遭受虫害(长足大竹象)时其相关研究尚未有报到.在本研究中，comp5294_seq2、comp7144_seq0、comp4915_seq1、comp1871_seq0、comp1486_seq1、comp6576_seq0、comp5390_seq0、comp26614_seq1、comp30333_seq1在长足大竹象取食之后表达量存在较大的变化.根据系统进化树的结果，comp1871_seq0、comp7144_seq0、comp6576_seq0与拟南芥WRKY70转录因子亲缘关系较近，而在拟南芥中的研究表明WRKY70是植物对病原菌响应过程中水杨酸(SA)和茉莉酸(JA)介导的信号事件的汇聚节点，且在决定依赖于SA和依赖于JA的防御途径之间的平衡中起着关键作用[47].对于comp1486_seq1，其与拟南芥的WRKY13的转录因子存在较高的同源性，而WRKY13在拟南芥中促进了茎的整体发育，因为WRKY13突变体的厚壁组织细胞数量、茎粗和维管束数量均减少.WRKY13突变体中木质素合成相关基因被抑制，导致茎秆较弱[48]，因此，comp1486_seq1在慈竹被长足大竹象取食时其表达量的上调可能与促进慈竹笋的发育来提升对虫害胁迫的抗性.除此之外，comp5390_seq0 与拟南芥的WRKY51转录因子及comp5294_seq2与拟南芥的WRKY60转录因子均存在较高的同源性, 而WRKY51蛋白介导了在低油酸(18:1)条件下JA介导的信号抑制[49]， WRKY60和WRKY40的共表达作用使植物对丁香假单胞菌和灰葡萄孢的抗性增强[50].根据上述结果，我们推测慈竹在遭受虫害的生物胁迫时，慈竹基因组内的WRKY转录因子可能参与了由SA合JA介导的信号转导(如comp1871_seq0、comp7144_seq0、comp6576_seq0、comp5390_seq0)及通过慈竹笋发育相关的WRKY转录因子(comp1486_seq1)表达量的上调来响应虫害的胁迫.

此外，本文从慈竹的虫害胁迫的转录组数据中鉴定了29个慈竹转录因子，其中26个WRKY转录因子具有完整的结构域，通过对其表达模式的分析可为以慈竹中WRKY转录因子如何响应长足大竹象取食胁迫及慈竹虫害的防控方面提供部分理论依据.