基于RNA-seq的福建柏R2R3-MYB转录因子鉴定和分析

周成城, 谢德金, 杨德明, 何天友, 陈凌艳, 郑郁善,

(1.福建农林大学园林学院；2.福建农林大学林学院, 福建 福州 350002)

转录因子(transcription factors, TFs)又称为反式作用因子，是一类能够与DNA序列上的启动子区域直接或间接地发生特异性作用的蛋白，从而激活或抑制目标基因的转录[1].转录因子家族参与植物生长发育、代谢等各个方面的基因表达调控，具有重要的生物学意义[2].MYB转录因子是植物中最早被克隆到的[3]，也是真核生物中功能最多样的家族[4].根据DNA结合结构域数目的不同，MYB转录因子分为1R-MYB、R2R3-MYB、3R-MYB和4R-MYB，其中R2R3-MYB是植物MYB转录因子家族中数量最多的一类[5].R2R3-MYB的主要功能是初级和次级代谢物的调节、植物生长发育调节、生物或非生物胁迫响应等[6].如拟南芥AtMYB83和AtMYB46是AtSND1的直接靶点，可诱导触发AtMYB58、AtMYB63和AtMYB85的表达，这些基因又通过与相关启动子AC元件相互作用来上调各种木质素合成基因[7-8]；R2R3-MYB的S7、S6和S2亚家族分别对黄酮醇[9]、花青素及原花青素[10]调控；在月季中，RhMYB108转录因子作为乙烯和茉莉酸信号接收器，通过激活和增强衰老相关基因表达从而调节花瓣衰老的启动的模式[11]；拟南芥AtMYB49通过调节角质层形成和抗氧化防御调节拟南芥的耐盐性[12].R2R3-MYB除了正调节作用，还可以作为负调控因子行使功能.在拟南芥中，AtMYB4能够抑制C4H的转录，从而抑制芥子酸酯的合成[13].以上研究均表明R2R3-MYB转录因子对植物的生长发育、物质合成等具有重要的作用.

福建柏[Fokieniahodginsii(Dunn) Henry et Thomas]为柏科福建柏属，为我国特有种属，分布在福建、广东、云南、广西等省区[14].其树干通直、树型优美、四季常青、适应性强、抗性较好、材性稳定、结构细腻，是我国优良的材用-景观两用型树种之一[15].木质素的合成对于木材的质量和树形有着重要影响，而R2R3-MYB转录因子对木质素合成的影响已经有了较为广泛的研究，如Du et al[16]对157个玉米MYB转录因子和125个拟南芥MYB转录因子进行分析，结果表明，有4个亚家族的MYB转录因子与木质素合成以及次生壁增厚有关，其中过表达AtMYB26能够使拟南芥木质部增厚；在菊花中，过表达CmMYB8使许多木质素合成组分编码基因表达下调，植物木质素含量降低[17]，同样地，银合欢的LlMYB1在转基因烟草中过表达也减少了木质素的含量[18].而目前福建柏在转录因子调控木质素合成的研究方面处于萌芽阶段，本文旨在利用福建柏转录组数据，对福建柏R2R3-MYB转录因子进行鉴定和生物信息学分析，为后续深入开展R2R3-MYB转录因子调控研究提供依据.

1 材料与方法

1.1 试验材料

试验材料为35年生生长良好、无病虫害的福建柏鳞叶、树皮(包括周皮和韧皮部)和根，于2019年8月采集于福建省福州市永泰县大湖国有林场的福建柏天然林，使用无菌水把材料清洗干净并用无菌滤纸擦干后，用锡箔纸包裹立即放进干冰中，带回实验室置于-80 ℃冰箱保存备用.3个部位材料在同一区域内选择年龄相近、长势相同的3个福建柏，各设置6个生物学重复.

1.2 试验方法

1.2.1 总RNA的提取和转录组测序 福建柏叶片、树皮和根的总RNA用总RNA试剂盒(天根生物有限公司)进行提取，提取后用琼脂糖凝胶电泳检测总RNA的质量和分光光度计(Nanodrop 2000)测定其浓度和纯度.并用反转录试剂盒(Takara)将总RNA反转录合成cDNA第一条链，置于-20 ℃冰箱中保存备用.从野外采集回来的样本立即送往武汉迈特维尔生物有限公司进行转录组测序，测序结果供后续分析使用.

1.2.2 福建柏MYB转录因子家族的筛选和鉴定 在福建柏转录组数据中，提取已有MYB基因家族功能注释的unigene数据，并提交至CD-HIT Suite在线数据库 (http://weizhong-lab. ucsd. edu/cdhit_suite/cgi-bin/index.cgi) ，去除冗余的核酸片段, 去冗余后的数据利用在线网站plantTFDB (http://planttfdb.cbi.pku.edu.cn/) 进行MYB转录因子筛选与鉴定.

1.2.3 福建柏MYB转录因子家族鉴定、理化性质与功能注释 利用OmicsBox分析工具对所有鉴定出的MYB转录因子序列进行氨基酸相似性比对、鉴定和筛选，并对筛选出的MYB转录因子进行GO功能分类.筛选出的MYB氨基酸序列N端所含DNA-Binding结构域使用在线预测工具ExPASy-pROSITE(https://prosite.expasy.org/)预测，并根据Dubos et al[4]的方法进行分类，分为1R-MYB(含MYB-related)、R2R3-MYB、3R-MYB、4R-MYB.运用在线处理软件ExPASy-Protparam (http://web.expasy.org/protparam/)、ExPASy-ProtScale (http://web.expasy.org/protscale/)和NetPho3.1 server (http://www.cbs.dtu.dk/services/NetPhos/)分析预测理化性质和亲疏水性.亚细胞定位预测使用在线网站Cell-PLoc(http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/)和WoLF PSORT(https://wolfpsort.hgc.jp/)，信号肽分析使用在线网站SignalP(http://www.cbs.dtu.dk/services/SignalP-4.1/)进行分析.

1.2.4 福建柏R2R3-MYB转录因子的保守结构域和进化树分析 利用MEME在线软件数据库(http://meme-suite.org/tools/meme)对福建柏R2R3-MYB的R2和R3保守基序进行分析.从在线网站plantTFDB数据库中下载拟南芥R2R3-MYB转录因子氨基酸序列，利用MEGA7.0软件与福建柏R2R3-MYB转录因子共同进行ClustalW序列比对和构建系统进化树.

1.2.5 福建柏R2R3-MYB转录因子的表达模式分析和验证 根据福建柏转录组数据中R2R3-MYB转录因子在叶片、树皮和根3个部位的FPKM(fragments per kilobase of transcript per million fragments mapped)值，利用R语言进行标准化和去中心化并使用heatmap包构建表达热图.选取不同表达模式下的R2R3-MYB转录因子进行RT-qPCR验证.利用Primer Premier 5.0软件进行引物设计(表1)，选择前期筛选出的Actin7作为内参基因，使用谢德金等[19]的方法进行qPCR扩增，每个样品做3个平行试验, 并做3个生物学重复.后用Ct(2-ΔΔCt)法计算福建柏R2R3-MYB转录因子在叶片、树皮、根中的相对表达量，利用Excel 2010和GraphPad 7.0进行统计、绘图.

表1 引物序列Table 1 Primer sequences

2 结果与分析

2.1 福建柏R2R3-MYB转录因子鉴定、理化性质与功能注释

在线网站plantTFDB初步筛选出145个福建柏MYB转录因子，经过氨基酸序列同源比对和保守结构域分析，有122个MYB转录因子具有MYB保守结构域且Blast同源性较高(均超过70%).通过分类得到71个R2R3-MYB、43个1R-MYB、6个3R-MYB和2个4RMYB，并依次对其进行编码FhMYB1～FhMYB122.71个R2R3-MYB的氨基酸数目为100～858个，相对分子质量11 618.52～96 926.56；总平均疏水系数(GRAVY)为-0.1～-1.0，推测都属于亲水蛋白；不稳定系数均高于40，都属于不稳定蛋白.Cell-PLoc和WoLF PSORT两个网站的预测结果显示71个R2R3-MYB转录因子均定位于细胞核中.信号肽分析表明全部序列的N端都没有信号肽剪切位点.

OmicsBox软件的GO分析结果显示,71个R2R3-MYB转录因子共注释到364个GO terms, 其中分子功能(molecular function, MF)、生物过程(biological process, BP) 和细胞组成(cellular component, CC) 涉及到的GO terms分别为150、82和132个.其中DNA binding类别中富集个数最多，在所有转录因子中均有注释，说明均含有DNA结合结构域；细胞组分中细胞结构组分(cellular anatomical entity)和胞内组分(intracellular component)均富集到66个R2R3-MYB转录因子，生物学过程中细胞学过程(cellular process)也相对较多，富集数目为50个.

2.2 福建柏R2R3-MYB转录因子家族保守结构域基序和进化树分析

R2-MYB和R3-MYB的特征保守基序为：-W-(X19)-W-(X19)-W-和-(F)-(X18)-W-(X18)-W-[4].由图1可知，在71个福建柏R2R3-MYB转录因子的R2-MYB特征中第13位、第33位和第53位氨基酸上出现色氨酸(W, Trp)，都分别间隔19个氨基酸，表明保守基序完全一致；在R3-MYB特征保守基序中，第13位苯丙氨酸 (F, Phe)被异亮氨酸 (I, Ile)取代，在第32位和第51位氨基酸上出现色氨酸，间隔18个氨基酸，符合R3-MYB保守基序特征.进一步分析表明，R2-MYB基序中第3个色氨酸残基(W, Tyr)前9位氨基酸残基的比对一致性高，表明该区域的氨基酸序列相对更保守. R3-MYB基序中2个色氨酸残基中间的18个氨基酸残基的比对一致性也较高，也表现出相对较高的保守性，如脯氨酸(P, Pro)、甘氨酸(G, Gla)、精氨酸(R, Arg)、苏氨酸(T, Thr)等.这些保守基序能够维持福建柏R2R3-MYB转录因子的螺旋—转角—螺旋(H-T-H)结构，进而与DNA更好地结合以行使功能.

R2：R2-MYB；R3：R3-MYB；W：色氨酸；F：苯丙氨酸；X：氨基酸.图1 福建柏R2R3-MYB转录因子的保守基序Fig.1 Conservative motif of the R2R3-MYB of F.hodginsii

基于生物学功能的研究，对拟南芥R2R3-MYB转录因子进行了系统的亚分类[4].将福建柏和拟南芥R2R3-MYB转录因子进行了多序列比对和构建系统进化树，并标注上拟南芥R2R3-MYB的亚类注释(图2).除亚类S12、S16、S18和S20外，其他亚类均有福建柏R2R3-MYB转录因子与拟南芥R2R3-MYB共同聚为一类，其中S5中聚类最多，有9个福建柏R2R3-MYB与AtMYB123聚为一类，可能参与苯丙烷类代谢途径的调控.有11个福建柏R2R3-MYB转录因子和拟南芥AtMYB111聚为一类，AtMYB111为S7分类，但根据拟南芥聚类方式，其未和AtMYB11和AtMYB12聚在一起，因此单独注释为W1，S7分类同样可能参与苯丙烷类代谢途径的调控.而有14个福建柏R2R3-MYB未和拟南芥R2R3-MYB聚为一类，而单独聚类为W2亚类.

图2 福建柏和拟南芥R2R3-MYB进化树Fig.2 Phylogenetic tree analysis of R2R3-MYB between F.hodginsii and Arabidopsis thaliana

2.3 福建柏R2R3-MYB转录因子基因的表达模式分析

对福建柏叶片、树皮和根的FPKM值进行表达模式热图分析(图3)，结果表明，福建柏R2R3-MYB转录因子的表达具有明显的组织特异性.根据表达模式的差异分为5个类别，其中分组Ⅰ包括FhMYB12、FhMYB23、FhMYB35、FhMYB36等16个转录因子，该类别在叶中表达与在树皮、根中差异显著，表达趋势为上调，且在树皮和叶的表达差异不大；分组Ⅱ包括FhMYB17、FhMYB19和FhMYB54，该类别的表达在根、树皮、叶中依次呈现阶梯式差异，表达趋势为上调；分组Ⅲ包括FhMYB10、FhMYB22、FhMYB31、FhMYB32等12个转录因子，该类别在树皮中表达与在叶、根中差异显著，表达趋势为下调，且在树皮和叶的表达差异不显著；分组Ⅳ包括FhMYB7、FhMYB8、FhMYB9、FhMYB11等21个转录因子，该类别在根中表达与在叶、树皮中差异显著，表达趋势为下调；分组Ⅴ包括FhMYB1、FhMYB2、FhMYB3、FhMYB4等19个转录因子，该类别在叶中表达与在树皮、根中差异显著，表达趋势为下调.

图3 福建柏R2R3-MYB表达模式Fig.3 Expression pattern of R2R3-MYB in F.hodginsii

为了进一步验证福建柏R2R3-MYB转录因子表达模式的准确性，从5种表达模式中选取具有代表性的转录因子(FhMYB8、FhMYB12、FhMYB17、FhMYB20和FhMYB38)，采用RT-qPCR的方法进行验证(图4)，结果表明，FhMYB8在根中相对表达量较高，在皮和叶中表达相对较低且两者差异较小；FhMYB12在根中相对表达量较高，在皮中次之，在叶中表达相对较低；FhMYB17在皮和叶中相对表达量较高，在根中相对表达量较低且与皮、叶中的表达差异较大；FhMYB20在叶中相对表达量较高，在根中接近无表达；FhMYB38在皮中相对表达量较高，3个部位的表达量差异相对较低.5个转录因子的RT-qPCR表达结果的总体趋势与FPKM值的表达模式基本一致，因此验证了R2R3-MYB转录因子的表达模式和组织特异性.

3 讨论与结论

R2R3-MYB转录因子是MYB转录因子家族中的一个大类，对植物的初生、次生代谢物的合成有着重要的作用[20].福建柏不仅材用价值、景观价值高，其精油的工业和药用价值也较高，这些特性都离不开福建柏初生、次生代谢物的合成，如木质素、黄酮类、萜类化合物等的合成.为探究福建柏中重要次生代谢物质的合成调控途径，对福建柏R2R3-MYB转录因子进行分析和研究是十分必要的.通过福建柏转录组数据库、plantTFDB数据库的筛选和比对，鉴定出71个福建柏R2R3-MYB转录因子，均含有R2、R3型MYB，且通过结构域分析和序列比对，R2-MYB、R3-MYB中基序中色氨酸残基前面的氨基酸残基的比对一致性高，表明该区域的氨基酸序列更加保守和集中.同时也说明 R2R3-MYB转录因子N端的保守重复序列能够更好地形成H(螺旋)-T(转角)-H(螺旋)结构,使R2R3-MYB转录因子蛋白与基因启动子区结合, 从而控制基因的转录水平[21].

在GO功能分类的分析中，预测福建柏R2R3-MYB转录因子主要富集于细胞过程和生物过程调控，这符合拟南芥[5]对R2R3-MYB转录因子的功能分类.因此进一步和拟南芥R2R3-MYB转录因子共同构建系统进化树，分析可知，除S12、S16、S18和S20外，其他亚类均有福建柏R2R3-MYB转录因子与拟南芥R2R3-MYB共同聚为一类，其中S5中聚类最多.在拟南芥中，S5亚家族参与花青素和黄酮类化合物合成途径[22].其次S13聚类中也较多，拟南芥S13中AtMYB61的研究表明，AtMYB61能够影响木质素的沉积[23].另外，W1类别中含有较多福建柏R2R3-MYB转录因子与AtMYB111聚在一类，AtMYB111原本在拟南芥中和AtMYB11、AtMYB12聚类，但在本研究中没有聚在S7中，推测可能由于构建进化树时，未将所有拟南芥R2R3-MYB转录因子进行聚类，而是将每个亚家族中的主要转录因子进行聚类.S7亚家族在拟南芥中有着非常重要的作用，不仅参与细胞氧化还原过程的调控[24]，也在黄酮类物质的调控中发挥作用[9].W2类别没有和任何拟南芥R2R3-MYB聚在一类的转录因子，该类别与S15亚家族在同一分支中，S15亚家族主要参与植物表皮细胞类型的分化调控[25-26]，W2类别是否符合S15亚家族的功能注释值得进一步探究.

基因在不同组织中的表达模式和基因的功能是密切相关的[27]，通过研究R2R3-MYB的表达模式，能够更好地理解R2R3-MYB在调控次生代谢物质合成过程的机理.FPKM表达模式分析表明福建柏R2R3-MYB转录因子具有明显的组织特异性，同时通过RT-qPCR进行了验证如：FhMYB8和FhMYB20分别在根和叶片中相对表达较高，且进化树分析都与AtMYB111聚为一类，AtMYB111属于S7亚家族，推测FhMYB8和FhMYB20可能分别在福建柏根和叶中参与黄酮类物质的合成；FhMYB38在树皮中表达较高，进化树分析与S15同一分支，推测FhMYB38可能在福建柏树皮中调控表皮细胞的分化.

综上所述，本研究基于福建柏叶、树皮和根的转录组测序数据，筛选和鉴定出71个福建柏R2R3-MYB转录因子.对71个R2R3-MYB转录因子进行了理化性质、GO富集分类、保守结构域、系统进化树分析以及表达模式的分析和验证，能够初步了解福建柏中R2R3-MYB的数量、功能分类以及表达趋势，对后续进一步开展福建柏R2R3-MYB转录因子对代谢物质的调控具有一定理论价值.