竹叶花椒MYB基因家族的鉴定及其表达特性的分析

冯发玉，王 毅，王丽娟，罗玛妮娅，徐令文，杨彦萍，陆 斌

（1.西南林业大学 林学院，云南 昆明 650224；2.云南省林业和草原科学院 a.云南省森林植物培育与开发利用重点实验室；b.国家林业局云南珍稀濒特森林植物保护和繁育重点实验室，云南 昆明 650201）

竹叶花椒Zanthoxylum armatumDC.为芸香科花椒属多刺小乔木。野生竹叶花椒主要分布在喜马拉雅山脉附近的亚热带和温带河谷地区，我国的四川、云南和重庆等西南地区已有广泛的栽植[1]。竹叶花椒的根、茎和果实均可入药，常用于祛风散寒、促进消化和行气止痛，其果实还是重要的香料和调味品，果皮花椒油的提取率高于红花椒，所含挥发油为杀虫剂和抑菌剂中的有效成分[2]。竹叶花椒具有早实丰产、抗病虫性强、耐干旱瘠薄、根系发达和固土保水保肥能力强等特点，是退耕还林工程中生态效益好的经济树种[3]。竹叶花椒是多刺的，其主干、枝、叶柄均密生皮刺，这可以保护竹叶花椒免受食草动物的攻击和机械伤害[4]，而其皮刺却不讨人喜爱，因为皮刺增加了栽培、管理、收获的难度[5]，提高了种植成本。因此，培育无刺的竹叶花椒品种是其育种计划的一个重要目标。然而，传统育种技术将无刺性状与其他重要性状结合起来，可能既耗时且所需费用又昂贵[6]。因此，挖掘竹叶花椒皮刺发育过程中的关键基因是十分重要的。

转录因子又称为反式作用因子，是由基因编码的一类调节蛋白，能够识别和结合特异的DNA序列，借助蛋白与DNA间的相互作用能有效调控目的基因在时空中的有序表达[7]。其中，MYB是植物体中最大的转录因子家族，对植物的生长发育具有重要调节作用，其位于N端的DNA结合区通常高度保守，由1～4个串联且不完全重复的R结构组成，且R结构是由50～53个氨基酸组成的[8]，位于羧基端的转录激活区，氨基酸一般折叠成双亲性，其一级结构的保守性不强，功能具有一定的可塑性[9]。有关研究结果表明，MYB类转录因子在一些表皮生长物（包括毛状体和皮刺）的发育中具有重要作用[6]。Zhao等[10]研究发现，拟南芥Arabidopsis thaliana毛状体的形成受到典型MBW转录复合物（其中包括MYB转录因子）的调控；Pandey等[11]在对热带刺茄Solanum viarum皮刺的研究中发现，与其皮刺发育有关的R2R3-MYB转录因子可能具有潜在作用；Swarnkar等[12]在对玫瑰Rosa rugosa皮刺的研究中发现，玫瑰树皮中表皮细胞的分化似乎是由MYB蛋白决定的，并且在其皮刺的发育过程中，参与木质素合成的MYB123和参与次生细胞壁生物发生的MYB43其表达水平均逐渐上调，其研究结果证实了成熟期间皮刺的木质化和硬化。

然而，到目前为止，尚未见有关于竹叶花椒MYB转录因子家族的研究报道。为了挖掘出与竹叶花椒皮刺发育有关联的MYB转录因子，研究其结构特点和生物学功能，本研究基于已有的竹叶花椒转录测序数据，对竹叶花椒中存在的MYB转录因子进行生物信息学及表达模式的分析，并结合已报道的其他物种中的功能性MYB转录因子，预测了部分竹叶花椒MYB转录因子的功能，筛选出可能参与皮刺发育的MYB转录因子，以期为进一步揭示MYB转录因子家族生物学功能提供理论依据，为后续无刺竹叶花椒品种的基因工程培育提供候选基因。

1 材料与方法

1.1 材 料

选取云南省林业和草原科学院树木园中生长状况良好的5年生实生竹叶花椒为研究材料。2019年4月单株取样，分别采集竹叶花椒的叶芽（样品编号为Ya）、根尖（样品编号为Gen）和皮刺刚突起（样品编号为PC1）、皮刺开始木质化（样品编号为PC2）、皮刺完全木质化（样品编号为PC3）等不同发育时期的皮刺及没有油点的果实样品（其编号为ZaF1）；2019年5月从相同样株上采集刚有油点的果实样品（其编号为ZaF2）。将所采样品迅速放入液氮中冷冻，送至派森诺生物科技股份有限公司，采用第二代测序技术（Next-generation sequencing，NGS），基于Illumina HiSeq测序平台，对7个样品进行双末端（Paired-end，PE）测序，测序读数为PE150，独立样本测序深度为6 G。

1.2 竹叶花椒MYB蛋白序列的获取

采用第二代测序技术获得竹叶花椒转录组数据，并从中筛选出含有MYB转录因子的数据，利用在线预测软件ORF Finder（http://www.ncbi.nlm.nih.gov/gorf/orfig.cgi）搜索获取具有全长氨基酸序列的竹叶花椒MYB转录因子序列。利用本地在线软件Blast比对分析MYB蛋白结构域，并剔除无MYB结构域的蛋白序列[13]，利用SMART在线软件筛选并确认竹叶花椒MYB蛋白序列是否包含SANT保守结构域。

1.3 竹叶花椒MYB蛋白特征分析

利用ExPASy protparam tool在线程序分析竹叶花椒MYB蛋白的理化性质，利用ProtScale在线程序分析蛋白的亲水性与疏水性，采用signaIP-5.0和PSORT Prediction在线软件预测蛋白的信号肽及亚细胞的定位，采用Prabi在线软件预测蛋白的二级结构[13]。利用MEME在线软件对39个竹叶花椒MYB蛋白的保守结构域进行预测，具体参数设置如下：基序位点分布情况，选择重复次数不限制；保守性基序的数目限制选择7，其他参数不变，均采用默认值[14]。

1.4 竹叶花椒MYB蛋白系统进化分析

以下载于Plant TFDB（http://planttfdb.gao-lab.org/）数据库的拟南芥MYB蛋白为参考序列，利用MEGA7软件构建竹叶花椒与拟南芥MYB转录因子蛋白系统进化树；将以不同组织间表达差异显著的8个竹叶花椒MYB转录因子蛋白序列与其他10个功能已被验证的植物MYB转录因子蛋白序列共同构建系统进化树，采用邻位连接法（Neighbor-joining，NJ）构建系统进化树，选用模型为Jones-Taylor-Thornton，执行方式设为Bootstrap method，自检重复次数设定为1 000次。

1.5 竹叶花椒MYB基因表达分析

以竹叶花椒的芽（Ya）、根尖（Gen）、不同发育时期的皮刺（PC1、PC2、PC3）和果实（ZaF1、ZaF2）的转录组测序（RNA-seq）数据为基础，获得竹叶花椒MYB转录因子基因表达谱数据，利用基因云在线软件上传基因表达谱数据，依据其FPKM值，在线绘制竹叶花椒MYB转录因子基因表达交互热图。

1.6 实时荧光定量PCR分析

为了验证转录组测序的结果及MYB转录因子的表达模式，设以特异引物TZaMYB10（F-5′GAATGTTCTCACAACATGAA3′、R-5′ATTAATAGTACTATTACTA3′）与TZaMYB26（F-5′GAGTTGCAGATTGCGATGGA3′、R-5′CCACGTCAGCAATCACAGGC3′） 和Ubiquitins（F-5′ CATGAAAGTCGAATATATCA3′、R-5′ ATGATTAACACTAGT3′）为内参基因。定量PCR反应体系为：2×SYBR的绿色主混合物12.5 μL，上、下游引物各0.5 μL，模板（cDNA）1 μL，ddH2O 10.5 μL，每个样本进行3次技术重复。使用PCR热循环仪（ABI 7300，应用生物系统，Foster City，CA，USA）对目标基因进行实时荧光定量检测分析，PCR的反应程序为：95 ℃下预变性5 min；95 ℃下变性30 s，58 ℃下退火30 s，72 ℃下延伸2 min（共35个循环）；72 ℃下终延伸7 min，4 ℃下保存。扩增完成后从65 ℃到95 ℃进行溶解曲线分析，验证扩增的特异性[15]。

1.7 ZaMYB10基因的克隆

以ZaMYB10基因序列为依据，通过引物设计软件Primer 5设计PCR特异扩增引物，引物由上海生工生物工程有限公司合成。克隆引物为ZaMYB10F（5′-ATGGGACGCCATTCTTCTTG-3′）和ZaMYB10R（5′-CATAAGTGGTGGAAGA-3′）。以竹叶花椒皮刺的cDNA为模板，分别以ZaMYB10F和ZaMYB10R为引物，以HiFiDNA聚合酶进行扩增，得到ZaMYB10全基因片段。经过电泳检测后，将目的片段连接到pUMT载体上，经PCR检测后送往上海生工生物工程有限公司测序。

2 结果与分析

2.1 ZaMYB蛋白的基本信息及理化性质分析

从竹叶花椒转录组数据中筛选出39条MYB转录因子序列，其所含氨基酸数量为101～621 aa，平均氨基酸数为288 aa。利用ExPASy protparam tool在线软件分析竹叶花椒MYB蛋白的基本性质，分析结果见表1。

由表1可知，在39个ZaMYB蛋白中，理论等电点（PI）的最大值为10.70，最小值为4.71，其均值为7.44，这一预测结果说明，竹叶花椒的MYB蛋白偏碱性；除ZaMYB6蛋白的亲水性的评分为正值外，其余38个ZaMYB蛋白的亲水性的评分均为负值，这一分析结果表明，ZaMYB蛋白是可溶性蛋白。39个MYB蛋白二级结构的预测结果表明，除了ZaMYB21、ZaMYB24和ZaMYB35蛋白的二级结构均以α-螺旋为主，其余36个ZaMYB蛋白的无规则卷曲的占比最高；总体而言，ZaMYB蛋白以无规则卷曲和α-螺旋为其主要结构，β-转角及延伸链均为其次要结构，无规则卷曲结构有利于ZaMYB蛋白稳定性的维持。信号肽的预测结果说明，竹叶花椒MYB蛋白不存在信号肽。亚细胞定位结果显示，39个ZaMYB蛋白位于细胞核的可能性最大，这一预测结果符合转录因子的分布特点。

2.2 竹叶花椒MYB转录因子保守结构的鉴定与分析

利用SMART在线程序对竹叶花椒MYB转录因子进行结构域分析，结果如图1所示。

图1 ZaMYB蛋白家族结构域的预测结果Fig.1 ZaMYB protein family domain prediction result

图1显示，39个竹叶花椒MYB蛋白具有共同的保守结构域，28个ZaMYB蛋白包含2个SANT结构域，其占比最高，这和大多数植物MYB转录因子结构域的分析结果相同，9个蛋白只包含1个SANT结构域，其中ZaMYB4和ZaMYB28除了具有1个SANT结构域外，均包含1个ZnF_C2H2锌指结构域，据此推测，其具有MYB转录因子功能之外的其他功能。利用MEME在线软件对竹叶花椒MYB家族蛋白系列进行分析，结果得到了如图2所示的7个保守元件，并命名其为基序1～7，此基序包含11～50个氨基酸，每个ZaMYB蛋白的基序数量为1～5个。在37个ZaMYB蛋白中发现了motif3，其为最常见的基序，其次分别是motif1和motif3，分别在32和30个ZaMYB蛋白中被发现。

图2 竹叶花椒MYB基因家族的保守元件和基因结构Fig.2 Conserved elements and gene structure of Z.armatum MYB gene family

2.3 竹叶花椒与拟南芥MYB蛋白的系统进化分析

拟南芥作为模式植物，有关其MYB转录因子基因功能的研究较为透彻，因此，结合拟南芥MYB蛋白的系统进化树进行分析，不仅有助于理解ZaMYB蛋白的系统进化关系，还可以推测ZaMYB蛋白的功能。本研究利用MEGA7软件构建竹叶花椒和拟南芥MYB转录因子蛋白的系统进化树，结果如图3所示。图3表明，竹叶花椒的ZaMYB36蛋白与拟南芥的AtMYB3R-1、AtMYB3R-2、AtMYB3R-3和AtMYB3R-4蛋白聚在同一分支中，均为R3类型；竹叶花椒的ZaMYB38蛋白可划分为R4类型，因为其与拟南芥的R4类蛋白聚在同一分支中。这一分析结果与SMART的结构域分析结果一致。整体而言，竹叶花椒和拟南芥的MYB蛋白主要聚为4个大分支：分支Ⅰ，包含16个R2R3-ZaMYB类蛋白和2个R1-ZaMYB蛋白；分支Ⅱ，包含4个R2R3-ZaMYB类蛋白和3个R1-ZaMYB蛋白；分支Ⅲ，包含8个R1-ZaMYB蛋白和1个R3-MYB蛋白；分支Ⅳ，包含4个无法与拟南芥MYB家族聚类的R3-MYB蛋白和1个R4-MYB蛋白，值得注意的是，在这4个无法聚类的蛋白中，ZaMYB4和ZaMYB28蛋白中均含锌指结构，而ZaMYB6蛋白具有未知结构域。依据拟南芥R1-MYB和R2-MYB蛋白的划分标准，可将ZaMYB蛋白进一步划分为23个亚组，其中，部分拟南芥亚组不包含ZaMYB转录因子。

图3 竹叶花椒和拟南芥MYB类蛋白的系统进化分析结果Fig.3 Phylogenetic analysis of MYB protein in Z.armatum and Arabidopsis thaliana

在拟南芥的蛋白中，S4和S7亚组成员主要参与黄酮类化合物和花青素的合成，因此推测，与其聚在一个分支的ZaMYB7、ZaMYB24、ZaMYB26和ZaMYB29蛋白可能都具有相似的功能[17]。根据竹叶花椒MYB转录因子基因表达谱数据，筛选出表达差异显著的8个ZaMYB转录因子，由这8个ZaMYB转录因子与拟南芥、甜橙Citrus sinensis和枇杷Eriobotrya japonica的功能已被验证的10个MYB转录因子共同构建的系统进化树如图4所示。由图4可知，竹叶花椒的ZaMYB10转录因子与拟南芥的AtMYB85转录因子和甜橙的CsMYB85转录因子均聚在一个分支中；竹叶花椒的ZaMYB25转录因子和拟南芥的AtMYB103转录因子聚在一个分支中，其bootstrap值为100。因此推测，竹叶花椒的ZaMYB10转录因子和拟南芥的AtMYB85、甜橙的CsMYB85转录因子可能都有类似的功能，竹叶花椒的ZaMYB25转录因子和拟南芥的AtMYB103转录因子可能具有类似的功能。

图4 部分竹叶花椒MYB蛋白和参与调控植物生长发育相关蛋白的进化关系Fig.4 Evolutionary relationship between MYB protein and proteins involved in plant growth and development of Z.armatum

2.4 竹叶花椒MYB转录因子的基因表达量分析

通过基因云在线软件上传竹叶花椒根、芽、皮刺和果的基因表达谱数据，依据其FPKM值，在线绘制如图5所示的竹叶花椒MYB家族转录因子基因表达交互热图。根据差异表达模式，将这些基因聚集成如下3个分支：分支Ⅰ，包含14个基因，除了ZaMYB5和ZaMYB33基因在果实中的表达量较高之外，其余12个基因在皮刺中的表达量均较高，其中的ZaMYB10和ZaMYB25基因在皮刺中均有特异表达；分支Ⅱ，包含11个基因，主要在根中表达；分支Ⅲ，包含14个基因，主要在芽和果实中有较高的表达量。值得注意的是，竹叶花椒MYB转录因子在皮刺不同发育时期的表达量逐渐升高，从图5中可以看出，其在皮刺发育起始时期（PC1）的表达量相对较低，而在皮刺发育木质化（PC3）时期的表达量较高，从皮刺刚突起到皮刺发育成熟阶段（从PC1到PC3），MYB转录因子的表达量总体呈现出逐渐增长的变化趋势，特别是ZaMYB10、ZaMYB23和ZaMYB25基因的表达量均呈现出明显的递增趋势，因此推测，这3个基因参与了对皮刺生长发育的调控。

图5 竹叶花椒MYB转录因子的基因表达情况Fig.5 Interactive heat map of gene expression of Z.armatum MYB transcription factor

2.5 竹叶花椒MYB基因的qRT-PCR分析

经特异引物检测反转录得到竹叶花椒cDNA，采用实时荧光定量PCR技术分析ZaMYB10和ZaMYB26基因在竹叶花椒的根、芽、皮刺、果实等组织中的具体表达情况，结果如图6所示。从图6中可以看出，ZaMYB26基因在发育早期的皮刺（PC1）和无油点果实（ZaF1）中的表达量均相对较高，ZaMYB10基因在皮刺中的表达特异，且其在发育早期（PC1）到成熟时期（PC3）的皮刺中的表达量呈现出逐渐增长的变化趋势，这一分析结果与转录组基因表达分析结果一致。

图6 ZaMYB10和ZaMYB26基因在竹叶花椒不同组织中的相对表达量Fig.6 Expression of ZaMYB10 and ZaMYB26 genes in different tissues of Z.armatum

2.6 ZaMYB10基因的克隆

以ZaMYB10F和ZaMYB10R为特异引物，通过RT-PCR反应扩增目的基因，PCR反应产物的琼脂糖凝胶电泳结果如图7所示。从图7中可以看出，EB显色后，在约为1 030 bp处有1条清晰的亮带，其与目的基因片段大小相符。切胶回收，将回收产物与pUMT载体连接，提取阳性克隆质粒，经BamHI和HindIII双酶切后得到的片段大小与PCR产物大小相吻合，利用NCBI ORF Finder软件对测序结果进行分析，结果表明，ZaMYB10基因具有全长为1 032 bp的完整的cDNA开放阅读框（ORF），其编码343个氨基酸。

图7 竹叶花椒ZaMYB10基因的克隆Fig.7 Cloning of ZaMYB10 Gene from Z.armatum

3 讨 论

在真核细胞中，基因表达调控是重要的生命活动，其中，转录调控主要是在一个复杂的由转录因子介导的相互作用系统中进行的。MYB转录因子普遍存在于植物中，对植物生长发育过程具有重要的调控作用[17]。随着相关研究的不断深入，MYB转录因子已在多个物种中被挖掘与鉴定，Hou等[18]在柑橘Citrus reticulata中发现了177个MYB转录因子；Dubos等[19]在拟南芥中发现了198个MYB转录因子；丁蒙蒙等[20]在楠木Phoebe zhennan中鉴定得到了82个MYB转录因子。本研究基于竹叶花椒转录组数据共挖掘出了39个MYB转录因子，其中包括9个R1-MYB、28个R2R3-MYB、1个R3-MYB和1个R4-MYB基因。对其理化性质的分析结果表明，不同的MYB转录因子编码的蛋白其理化性质存在差异，说明竹叶花椒不同的MYB转录因子发挥的生物学作用也不相同[21]。亚细胞定位预测结果显示，39个MYB转录因子定位于细胞核的可能性最大，这一预测结果符合转录因子调控下游基因表达的特性[22]。

竹叶花椒主干和枝条均密生皮刺，对果实采摘造成了一定的困难，故采摘果实常需大量的劳动力，费工费时[5]。目前，关于竹叶花椒皮刺的研究相对薄弱，因此，通过分析3个不同发育时期的皮刺中MYB转录因子的表达情况，以挖掘出可能参与皮刺生长发育的MYB转录因子。有关研究结果表明，在植物不同发育时期的同一组织中，基因表达存在差异的MYB转录因子在其发育过程中起着重要作用。Jia等[23]在对甜橙的研究中发现，CsMYB308和CsMYB330基因的表达水平在甜橙果汁囊的木质化过程中发生了显著变化，他们通过瞬时分析和过量表达分析证明了这两个基因在甜橙果汁囊的木质化中均起着重要作用；李志根等[24]在分析MaMYB2基因在不同组织中的表达模式后发现，MaMYB2基因可能参与了葡萄风信子Muscari botryoides花青素积累的转录调控。研究中发现，ZaMYB10、ZaMYB23和ZaMYB25基因的表达量从早期皮刺（PC1）到成熟皮刺（PC3）逐渐增长，据此推测，这3个基因可能参与了皮刺生长发育的调控。

Asano等[25]在对玫瑰皮刺形态的研究中发现，木质素的积累从皮刺顶部开始，向下扩散到皮刺的整个部分中，这表明皮刺的硬化与木质素的积累密切相关。竹叶花椒的皮刺和玫瑰的皮刺具有相同的发育史，其皮刺形态均经历了从小到大和从软到硬的变化过程。最近有关研究者[12]发现，玫瑰皮刺的发生和毛状体的形成类似，受到典型MBW转录复合物的调控，其中皮刺的形成似乎是由MYB蛋白决定的，玫瑰皮刺在其形成过程中积累了大量的次生代谢产物，尤其是木质素和类黄酮物质，一些与木质素合成相关的MYB转录因子，在皮刺发育过程中的表达水平逐渐上调，这一结果表明，MYB转录因子在皮刺的发育过程中具有重要作用。在对竹叶花椒MYB转录因子系统进化关系的分析中发现，竹叶花椒的ZaMYB10基因与拟南芥的AtMYB85基因和甜橙的CsMYB85基因均聚在一个分支中，表明竹叶花椒的ZaMYB10与拟南芥的AtMYB85和甜橙的CsMYB85的亲缘关系均较近，其ZaMYB25和拟南芥的AtMYB103也聚在一个分支中，其bootstrap值为100，且拟南芥的AtMYB85[26]、AtMYB103[27]基因和甜橙的CsMYB85[28]基因均参与了木质素的生物合成。结合基因表达分析结果推测，ZaMYB10和ZaMYB25基因在皮刺发育过程中均发挥着重要作用。

为了验证测序数据的准确性，采用实时荧光定量PCR分析发现，ZaMYB10基因在皮刺中的特异表达情况和转录组分析结果一致，由此推测，ZaMYB10基因可能参与了竹叶花椒皮刺中木质素的生物合成，与皮刺的硬化具有密切的关系。除此之外，对ZaMYB10基因的克隆结果表明，ZaMYB10基因具有全长为1 032 bp的完整的cDNA开放阅读框（ORF），其编码343个氨基酸，这进一步证实了该基因的真实性。这一研究结果为进一步探讨竹叶花椒MYB转录因子参与皮刺生长发育调控的机理奠定了基础，同时为无刺或少刺的竹叶花椒的基因工程培育提供了候选基因。然而，本研究对ZaMYB10基因的功能未加验证，后续将通过异源表达和基因敲除等验证其功能。

4 结 论

从竹叶花椒中共鉴定出39条MYB转录因子序列，其均包含SANT保守结构域，其中包括9个R1-MYB型基因、28个R2R3-MYB型基因、1个R3型基因和1个R4型基因。分析其系统进化关系和基因表达情况后发现，ZaMYB10和ZaMYB25基因可能与其皮刺的发育有关系，且实时荧光定量PCR分析结果也表明，ZaMYB10基因在皮刺中有特异表达，其表达模式和转录组分析结果一致；此外，本研究成功克隆了ZaMYB10基因。这些证据证明，ZaMYB10基因真实存在并且可能在皮刺发育过程中具有重要作用。