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

张鑫月，蒋筱蝶，李沅秋，罗朝兵，，唐 昊，

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

转录因子（Transcription factor，TF）通过目标基因的启动子区域的结合在激活基因表达的过程中发挥着至关重要的作用。基本亮氨酸拉链（bZIP）结构域，包含一个碱性氨基酸区域和一个由60～80个氨基酸残基构成的亮氨酸拉链[1-2]，是真核生物中最大的转录因子家族之一，通常与植物的生物或非生物胁迫相关，参与了植物中许多重要的生理生化过程[3]。碱性氨基酸区域有一个保守的N-x7-R/K基序，负责DNA的结合和核定位，而亮氨酸拉链区由多个亮氨酸残基及其它疏水性氨基酸（如异亮氨酸、缬氨酸、苯丙氨酸或甲硫氨酸等）形成的7次重复[4]。

目前，通过使用全基因组分析，在多种植物中预测并确定了大量的bZIP转录因子，如在拟南芥[1]中发现了75个bZIP转录因子，在水稻[5]中鉴定的89个bZIP转录因子，在高粱[6]中鉴定的92个bZIP转录因子，玉米[7]中鉴定的125个bZIP转录因子，二穗短柄草[8]中鉴定的96个bZIP转录因子，小麦[9]中鉴定的187个bZIP转录因子等。此外，大量的研究表明，bZIP转录因子在许多组织和器官的发育方面具有重要作用，如维管束的发育、花的诱导和发育[10-12]、种子的成熟和萌发[13-14]。同时，其他研究也报道了bZIP转录因子在脱落酸（ABA）[15]、赤霉素（GA）[16]、乙烯[17]激素水平引起的调控和对病原体感染[18]、干旱[19]、寒冷[20]、高温[21]和高盐[22]的应激反应中起关键作用。

慈竹（Bambusa emeiensis）是丛生竹类，广泛分布于我国南方地区，被广泛应用于造纸及园林绿化等领域[23]，如在造纸领域，按照竹种对造纸质量的影响进行分级排序，竹种可以被分为4级，而慈竹为第一级[24]。长足大竹象（Cyrtotrachelus buquet）属鞘翅目象甲科，是主要的竹林害虫之一。成虫期的长足大竹象可将喙完全伸入幼嫩竹笋内部，取食竹腔内层[25]，引起竹腔变黑、腐烂，造成新生竹发育不良，形成脓包竹及造成大量退笋，这将极大地影响了慈竹的产业化发展[26-27]。

目前尚未有慈竹基因组草图公布和慈竹在虫害（长足大竹象）生物胁迫下bZIP转录因子的分析，因此结合虫害（长足大竹象）取食胁迫下慈竹转录组的测序结果及生物信息学分析，对慈竹的bZIP转录因子进行鉴定及分析，以期发现bZIP转录因子在慈竹抗虫害表征中所发挥的作用，同时为研究慈竹bZIP转录因子的结构和功能提供参考。

1 材料与方法

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

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

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

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

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

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

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

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

在Plant Transcription Factor Database（http://planttfdb.cbi.pku.edu.cn/.）中下载山羊草（Aegilops tauschiia）bZIP蛋白序列，使用DNAMAN6.0及MEGA 7.0软件对慈竹及山羊草的bZIP蛋白进行多序列比对，然后在MEGA 7.0软件中的系统进化树构建方法中选用邻接法并设置1000次重复bootstrap检验、Poisson模型构建bZIP蛋白的系统发生树[32-33]。利用在线网站（http:// meme-suite.org/index.html）并使用默认参数设置[34]对bZIP蛋白进行motif的分析。

1.5 慈竹bZIP表达模式

根据慈竹中鉴定出的bZIP转录因子的相对表达丰度值，将值导入MeV 4.9.0分析软件中，制作热图并进行聚类分析[35]。

2 结果分析

2.1 bZIP蛋白的性质分析

根据ORFFinder查找结果，筛选出10个大于100个氨基酸的ORF序列。这些bZIP蛋白的ORF的氨基酸含量差别较大，最长ORF的氨基酸数从149（comp719_seq1）到641（comp9572_seq1）；相对分子量位于16221.19 Da （comp719_seq1）到68353.04 Da（comp9572_seq1）之间；等电点介于 4.56（comp4125_seq0）～9.73（comp10711_seq0）之间；此外，蛋白质的亲水性总平均值结果表明慈竹bZIP蛋白均为亲水性蛋白质（表1）。

表1 慈竹bZIP蛋白的相关性质

图1 慈竹及部分山羊草bZIP转录因子的序列比对图

2.2 慈竹bZIP蛋白的结构域分析

bZIP结构域包含有2部分，一个是包括核定位信号以及可以结合DNA序列的N-x7-R/K基序在内的一段碱性区域，另一段是亮氨酸拉链结构域，具有L-x6-L-x6-L结构，参与bZIP蛋白与DNA的二聚化[4]。从NCBI下载6条已知的慈竹bZIP转录因子氨基酸序列，并从Plant Transcription Factor Database中挑选出20条山羊草bZIP转录因子氨基酸序列，利用MEGA7.0和DNAMAN软件对这36条bZIP转录因子进行多序列比对。在这10个新鉴定的慈竹bZIP转录因子氨基酸序列中，它们的碱性氨基酸区域内的N-x7-R/K位点具有极高的保守型，说明该结构域在bZIP在结合启动子区上起关键的作用。10个新鉴定的bZIP转录因子中，有9个bZIP转录因子具有完整的结构域，而comp10711_seq0在核定位信号存在一定的残缺。此外，对bZIP转录因子的氨基酸残基进行分析后发现，在碱性氨基酸区后接的7重重复的亮氨酸拉链区域，前3个重复的亮氨酸（L）在相应位点上较保守，但在少数情况下会被缬氨酸（V）和甲硫氨酸（M）等替代（图1）。

2.3 慈竹及山羊草bZIP转录因子的系统进化树构建

从NCBI的protein数据库中以“Bambusa emeiensis bZIP”为索引下载的6条已知的慈竹bZIP序列，20条从Plant Transcription Factor Database中下载的山羊草bZIP转录因子序列及10条新鉴定的慈竹bZIP转录因子，利用MEGA 7.0软件中CLUSTAL W功能对这36条bZIP蛋白序列进行比对，绘制出慈竹bZIP蛋白的进化树（图2）。根据系统进化树的分析结果，可以将36条序列分为3组。Group 1含有19条序列（包括慈竹comp10711_seq0、BebZIP4、comp4125_seq0等8条慈竹bZIP转录因子），Group 2含有包括comp8646_seq0、BebZIP3在内的4条序列。Group 3共有13条序列（包含BebZIP2、comp_seq1和comp2588等在内的6个慈竹bZIP转录因子）。

图2 慈竹及山羊草bZIP转录因子的系统进化树

2.4 慈竹bZIP转录因子motif分析

利用MEME软件分析慈竹bZIP转录因子的motif（图3），结果发现：有15个慈竹bZIP转录因子含有motif1，结合motif1的氨基酸代表性序列，可推测其主要包含核定位信号，可以结合DNA序列的N-x7-R/K基序在内的一段碱性区域及部分亮氨酸拉链结构域，这与多序列比对结果相吻合（图1）。此外，有1个慈竹bZIP转录因子（comp10711_seq0）未检测到motif1，但其与已经鉴定的慈竹BebZIP4的序列对比，可知comp10711_seq0 bZIP结构域存在残缺，可能是由于转录组数据拼接不完整导致的序列不完整。

图3 慈竹bZIP转录因子motif分析

2.5 慈竹bZIP基因表达分析

根据本试验的结果（图4）：在不同竹笋的未被取食部位（Group1与Group2），bZIP转录因子在慈竹中表达模式基本相同，而在同一根竹笋的取食与未被取食部位（Group1与Group3），而在经长足大竹象取食前后出现表达量的明显变化，初步说明bZIP转录因子可能参与了慈竹对虫害的响应。其中comp8738_seq0、comp10711_seq0、comp4125_seq0、comp2588_seq0，在未取食部位与取食部位，表达量变化较为明显，说明这几个bZIP转录因子可能对慈竹抵御虫害有重要作用；此外，comp719_seq1在3组中含量均高，表明可能与竹笋的生长发育有关。

图4 慈竹bzip基因的表达模式

3 讨论与结论

bZIP转录因子在植物响应生物与非生物胁迫的应激反应中具有非常重要的作用[18-22]。通过生物信息学手段对慈竹的转录组数据进行分析，从中新鉴定出10个慈竹bZIP转录因子。

通过系统进化树上的聚类分析，发现BebZIP3和comp8646_seq0，BebZIP4和comp10711_seq0同源性较高；同源性较高的基因，其功能也可能具有相似性。在早前的慈竹bZIP研究中已鉴定出6个慈竹bZIP转录因子，其中BebZIP2、BebZIP5和BebZIP6转录因子参与了盐、干旱和ABA引起的非生物胁迫应激响应[36-38]，BebZIP1转录因子可能参与了慈竹中组蛋白基因的表达，进而影响其功能，从而改变慈竹的生长发育的调控[36，39-40]；BebZIP3转录因子可能有助于提高慈竹对多种非生物胁迫的适应能力[36，41]，而BebZIP4的调控功能尚不清晰。此外，根据长足大竹象取食慈竹笋前后转录组数据结果，10个新鉴定的慈竹bZIP转录因子在取食处理前后，表达量各有差异。其中comp8738_seq0、comp10711_seq0、comp4125_seq0、comp2588_seq0等转录因子表达量差异极为明显，推测这几个慈竹bZIP转录因子与长足大足象取食导致的慈竹应激胁迫响应相关。comp2588_seq0在受到长足大竹象取食胁迫后表达量下调，其具体原因还有待进一步分析。comp8738_seq0、comp10711_seq0、comp4125_seq0表达量上调，且它们和山羊草的EMT13093、EMT18702具有较高的同源性，而与这2个山羊草同源的水稻bZIP转录因子参与了内质网的应激反应，同时作为内质网压力传感器，激活转录因子BZIP50和伴侣BIP1[42-43]，因此推测comp8738_seq0、comp10711_seq0、comp4125_seq0可能也具有类似的调控作用，即通过调节细胞内内质网的应激反应来应对长足大竹象的取食胁迫。

本研究通过对慈竹bZIP转录因子的系统进化树分析、蛋白质保守结构域的预测及在长足大竹象取食胁迫下bZIP转录因子表达量变化的分析，将为慈竹bZIP转录因子在响应虫害方面的研究提供一定的借鉴意义。