玉米MYB转录因子靶基因的全基因组预测及验证

张晓芳，于好强，罗 羲，张于福，李晚忱，付凤玲

(四川农业大学玉米研究所/农业部西南玉米生物学与遗传育种重点实验室,四川 成都 611130)

转录因子对靶基因启动子结合位点的识别与结合，是基因表达调控最关键的环节，也是功能基因组学分析的重要内容[1,2]。以往通过酵母单杂交、染色质免疫共沉淀(ChIP)、逆转录PCR(RT-PCR)、基因表达芯片等实验手段筛选转录因子特异识别的启动子，鉴别其靶基因，准确性不高，效率低。随着生物信息学数据的积累，结合机器学习等智能计算技术的应用，根据转录因子结合位点以及转录因子与其靶基因表达模式的相关性，已开发出一些预测转录因子靶基因的生物信息学软件[3-5]。

MYB是人类成髓细胞瘤(MYB)转录因子家族的植物直系同源蛋白，根据其包含MYB重复结构的个数分为不同亚家族。在植物中，大部分MYB转录因子均包含R2和R3两个MYB重复结构，因此划分为R2R3-MYB亚家族，调控植物众多生长发育及逆境应答相关基因的表达[6-12]。玉米中已克隆的R2R3-MYB亚家族基因有C1、P1、MYB-IF25、MYB-IF35等，都与逆境应答有关[11-12]。但是，对于该亚家族更多MYB转录因子的靶基因及调控机制，却鲜有报道[13]。R2R3-MYB亚家族转录因子识别启动子的核心序列为TAACTG，其中第三个碱基A具有高度保守性，在MYB识别靶基因中起关键作用。但是，启动子具有此核心序列的基因并不全是R2R3-MYB转录因子调控的靶基因。

为了鉴定更多R2R3-MYB转录因子调控的靶基因，进一步研究其在玉米逆境应答及生长发育过程中的功能，本研究根据MYB识别序列的核心序列及其侧翼序列，结合运用HexDIFF算法和支持向量机(SVM)构建分类模型，在玉米全基因组范围内对MYB转录因子的靶基因进行预测，根据功能注释分析其可能的生物学功能。以电泳迁移率实验(EMSA)在体外验证MYB转录因子与预测序列的结合，并在玉米愈伤组织中瞬时表达检测预测的MYB靶基因启动子的活性。

1 材料和方法

1.1 模体框数据集构建

综合运用加州大学圣克鲁兹分校开发的基因组浏览器UCSC Genome Browser(http://genome.ucsc.edu/)、真核启动子数据库EPD(http://www.epd.isb-sib.ch/)、转录因子数据库TFD(http://www.ifti.org/)、转录调控区数据库TRRD(http://wwwmgs.bionet.nsc.ru/mgs/gnw/trrd)、真核转录因子基因组结合位点和DNA结合谱数据库TransFac(http://www.gene-regulation.com/)、转录因子与DNA结合位点模体数据库Jaspar(http://jaspar.cgb.ki.se/cgi-bin/jaspar_db.pl)、转录起始位点数据库dbTSS(http://dbtss.hgc.jp/)、植物顺式作用调控元件及启动子序列分析数据库PlantCare(http://oberon.fvms.ugent.be:8080/PlantCARE/index.html)、植物顺式调控DNA元件中的基序数据库PLACE(http://www.dna.affrc.go.jp/htdocs/PLACE/)、植物启动子数据库PlantPromDB(http://mendel.cs.rhul.ac.uk/mendel.php?topic=plantprom)、美国国立生物技术信息中心文献数据库PubMed(www.ncbi.nlm.nih.gov/pubmed)等数据库，搜索经酵母单杂交、ChIP、RT-PCR、基因表达芯片等技术筛选鉴定的植物R2R3-MYB转录因子的靶基因。用PERL语言编写脚本，分别从拟南芥资源网(www.arabidopsis.com)、禾本科资源网(www.gramene.org)和玉米基因组数据库(www.maizesequence.org)下载转录起始位点上游2000 bp序列，用启动子预测软件Promoter 2.0(http://www.cbs.dtu.dk/services/Promoter/)预测其中的启动子序列，分别构建Arab+orysa_promoter和Zea_promoter数据库，从中截取已验证MYB转录因子结合位点核心序列TAACTG及其两侧各100 bp的序列，构建阳性模体框数据集(Positive motif frame data)。从Arab+orysa_promoter数据库中，检索已验证MYB转录因子以外但具有核心序列TAACTG的序列，截取其两侧共206 bp序列，构建阴性模体框数据集(Negative motif frame data)；以TAACTG核心六聚体为中心，在两侧以A、T、C、G四种核苷酸为元素随机函数生成100 bp的序列，构建随机模体框数据集(Random motif frame data)；再从Zea_promoter数据库中，检索核心序列TAACTG，截取其两侧共206 bp序列，构建待测模体框数据集(Test motif frame data)。前3个数据集又合称为训练模体框数据集(Train motif frame data)。

1.2 玉米MYB转录因子识别启动子预测

用PERL语言编写脚本，在上述各模体框数据集序列中截取所有可能的6 bp六聚体，用HexDIFF算法分别计算其在阳性模体框数据集、阴性模体框数据集和随机模体框数据集中出现的频率R(h)[14]：

参照Huang等介绍的方法[16]，将阳性模体框数据集、阴性模体框数据集和随机模体框数据集中所有可能六聚体的出现频率，用智能计算的SVM，综合使用线性和非线性方法分类。将训练模体框数据集和待测模体框数据集归一化处理，在一定的保证概率预测待测模体框数据集中R2R3-MYB转录因子识别的核心序列TAACTG及侧翼序列，用WEBLOGO软件(http://weblogo.berkeley.edu/logo.cgi)在线生成模体框碱基频率图[17]。

1.3 基因功能富集分析

用禾本科资源网(http://www.gramene.org)的GmeneMART在线软件，选择Ensembl Gene ID滤器，从Plant Genes的Zeamaysgenes(AGPv2)数据库，导出所预测MYB转录因子靶基因的功能注释。再将导出的靶基因功能注释(GO)导入AgriGO在线软件的Singular Enrichment Analysis(SEA)工具[18]，选择玉米基因组ZeamaysAGPv2作背景，以0.05显著水平、Fisher统计测试算法提交数据，得到靶基因的功能富集分析结果。

1.4 凝胶迁移率检测

从玉米模式自交系B73的cDNA 样品扩增R2R3-MYB转录因子基因ZmMYB-IF25的开放阅读框序列，经DNAworks在线软件(https://omictools.com/dnaworks-tool)对密码子偏好性优化后，构建原核表达载体pET-21a(+)-ZmMYB-IF25，转化大肠杆菌BL21菌株，在0.4 mmol/L异丙基-β-D-硫代半乳糖苷(Isopropyl-β-D-thiogalactoside, IPTG)诱导下，37℃ 240 r/min 培养3 h，用12%十二烷基硫酸钠聚丙烯酰胺凝胶电泳(Sodium dodecyl sulfate polyacrylamide gel electrophoresis, SDS-PAGE)进行表达检测后，装柱洗脱，分离纯化MYB-IF25蛋白。

随机选取30条预测模体框序列(表1)，生工生物工程(上海)股份有限公司合成，以已验证的玉米靶基因Bz1和Bz2启动子MYB核心序列TAACTG及两侧各106 bp序列为阳性对照，用地高辛标记后与纯化的ZmMYB-IF25蛋白在5×EMSA结合缓冲液中混合，以2倍浓度的未标记双链探针为冷竞争结合对照，室温保持20 min上样，进行连续非变形聚丙烯酰胺凝胶电泳，转膜和紫外交联后化学发光并成像，根据显色条带的迁移率鉴别MYB-IF25蛋白与预测模体框序列的结合。

1.5 瞬时表达验证

从体外验证的预测模体框序列中随机选取4条，扩增其全长启动子序列，构建成启动β-葡萄糖苷酸酶(β-glucuronidase)基因GUS的瞬时表达载体pBI221-pCandidate-GUS，混合1/3荧光素酶(Luciferase)基因LUC表达载体pUbi-LUC为内参，用DuPont PDS1000/He型基因枪(Bio-Rad, USA)，在1100 psi系统压力和25 mm汞柱真空度下，转化玉米胚性愈伤组织。每个处理重复3次，每次重复转化15块大小、形状和质地基本一致的愈伤组织。转化后在8%甘露醇高渗培养基上27℃黑暗培养24 h，GUS组织化学染色后在解剖镜下观察拍照蓝色斑点。剩余愈伤组织用荧光素酶检测试剂盒(Promega，上海)显色后，用Fluoroskan Ascent FL型荧光/化学发光微孔检测仪(天齐生物科技有限公司，上海)，检测高渗胁迫0(对照)和24 h的GUS与LUC荧光强度比值GUS/LUC，对预测靶序列响应高渗胁迫的体内启动活性进行相对定量。

表1 随机选取的预测模体框序列

2 结果与分析

2.1 玉米MYB靶基因及启动子结合位点

当SVM参数C=2，γ=0.03125时，以98.38%的保证概率，从Zea_promoter数据库136770条序列中，预测到435个MYB转录因子结合位点，对应下游424个靶基因。其中，Bz1(GRMZM2G165390)和Bz2(GRMZM2G016241)是已验证的MYB转录因子靶基因[19]。与玉米基因组序列比对结果，绝大部分预测靶基因随机分布于全部10个玉米染色体上，没有染色体偏好，但有的区段有集中分部情况，没有分布于线粒体和叶绿体基因组，极少部分预测靶基因不能与基因组序列匹配(图1)。

预测靶基因启动子MYB识别位点及其侧翼序列的碱基频率与已知MYB靶基因相似，具有高度保守性，但距离核心序列较远的碱基存在一定的变异(图2)。这既与物种差异有关，也可能是同一物种基因间的差异。

图1 预测MYB靶基因在玉米染色体上的分布Fig.1 Distribution of predicted MYB targeted genes on the maize chromosomes

2.2 玉米MYB靶基因的生物学功能

用GrameneMART软件从Plant Genes的Zeamaysgenes(AGPv2)数据库中共搜索到229个预测MYB靶基因的功能注释共计1471条。其中，45个(19.7%)基因涉及刺激响应(GO：0050896)，33个(14.4%)基因介导逆境应答(GO：0006950)，10个(4.3%)基因与再生过程(GO：0000003)相关，10个(4.3%)基因参与发育过程(GO：0032502)。以细胞组分类分析发现，93个(40.6%)基因编码蛋白定位于细胞内(GO：0005622)，12个(5.2%)基因编码蛋白定位于细胞外(GO：0005576)。

2.3 MYB转录因子与预测靶基因启动子序列的体外结合

EMSA结果表明，随机选取的30条预测模体框序列中，除基因ID为GRMZM2G086773、GRMZM2G137596和GRMZM2G087719的3个基因的预测模体框序列以外(表1)，其余27条(90%)预测模体框序列均可在体外与MYB-IF25蛋白结合，电泳迁移条带滞后(图3)，表明本研究所用的转录因子全基因组预测方法具有较高的可靠性。

图2 MYB靶基因核心及侧翼序列碱基频率Fig.2 Base frequency during core and flanking sequences of MYB targeted genes

FP: 无MYB-IF25蛋白结合的自由探针；B：MYB-IF25蛋白与地高辛标记预测模体框序列结合；C：MYB-IF25蛋白与未标记预测模体框序列结合(冷竞争结合)；1～30：随机选取的30条预测模体框序列；bz1和bz2：Bz1和Bz2基因核心序列TAACTG及两侧106 bp序列(阳性对照)。FP: free probes without the MYB-IF25 protein combination; B: the predicted motif frame sequences marked with gidoxin and combining to the MYB-IF25 protein; C: the predicted motif frame sequences unmarked but combining to the MYB-IF25 protein (cold target competition); 1～30: randomly selected samples of the predicted motif frame sequences; bz1 and bz2: the core sequence and its 106 bp flanking sequences of the Bz1 and Bz2 genes (positive control).图3 ZmMYB-IF25蛋白与预测模体框序列结合的EMSA检测Fig.3 EMSA detection for combination between the ZmMYB-IF25 protein and the predicted motif frame sequences

2.4 预测靶基因启动子序列的体内表达活性

从体外验证的27条预测模体框序列中随机选取4条涉及逆境刺激响应相关靶基因(GRMZM2G044829、GRMZM2G115698、GRMZM2G310161和GRMZM2G430675)，扩增其全长启动子序列，构建成启动GUS基因的瞬时表达载体pBI221-pCandidate-GUS，基因F枪法转化玉米愈伤组织，8%甘露醇高渗培养基27℃黑暗培养24 h后，GUS组织化学染色结果显示，4个预测靶基因的启动子均能在玉米愈伤组织内启动GUS基因的表达，产生显色反应(图4)。这4个启动子8%甘露醇渗透胁迫下启动GUS表达产生的GUS荧光强度与LUC荧光强度的比值，均比胁迫前显著增高(图5)，说明这4个靶基因的启动子在玉米细胞内具有启动活性，且受高渗胁迫诱导。

3 结论与讨论

本研究将HexDIFF算法与SVM结合，用已知MYB结合位点序列建模，在玉米全基因组范围内预测到424个MYB靶基因和435个MYB结合位点(图1、2)，涉及逆境刺激响应、再生和发育等众多生长发育过程。EMSA实验表明，预测的结合位点与玉米MYB可相互结合(图3)。GUS瞬时表达证实预测得到的MYB靶基因启动子具有启动活性(图4、5)，而且预测的逆境相关靶基因启动子启动报告基因GUS，在渗透胁迫下的表达量显著提高，更说明预测的MYB靶基因有一部分的功能确与植物抗逆反应相关。以上结果证明本研究预测方法可靠性高，为预测转录因子靶基因提供借鉴。随着第三代测序技术(尤其是转录组学)的快速发展，人们必将对玉米MYB转录因子进行更为深入的研究，届时，玉米MYB转录因子参与的生长、发育等调控机制将会被更好地解析。

8：GRMZM2G044829，19：GRMZM2G115698，29：GRMZM2G310161，30：GRMZM2G430675图4 MYB靶基因启动子启动GUS基因瞬时表达Fig.4 Transient expression of the GUS gene under the control of the promoters of the predicted MYB target genes

8：GRMZM2G044829，19：GRMZM2G115698，29：GRMZM2G310161，30：GRMZM2G430675图5 MYB靶基因启动子在高渗胁迫下的启动活性Fig.5 Promotion activity of the promoters of the predicted MYB targeted genes under osmotic stress

EMSA体外验证实验发现，预测结果可能存在10%假阳性，可能的原因是各数据库能够检索到的MYB靶基因及其结合位点有限，构建的阳性模体框数据集代表性可能不够全面。随着转录因子靶基因数据量的积累，本研究方法可能会得到更好的验证。本研究的分析还表明，核心序列TAACTG是MYB结合位点的必要条件，在我们检索的阳性数据和预测的全部玉米MYB结合位点中100%保守，但并不是充分条件，MYB的的结合还决定于核心序列的侧翼序列(图2)。用TRANSFAC(http://gene-regulation.com/pub/databases.html)和JASPAR(http://jaspar.genereg.net/)等软件搜索可发现， TAACTG序列还存在于DREB(Dehydration responsive element binding)、NAC(NAM/ATAF1/2/CUC2)等转录因子的启动子序列中[20-21]。这也说明，MYB对其靶基因的转录调控可能还与DREB、NAC等转录因子存在协同关系。

参考文献：

[1] Levine M, Tjian R. Transcription regulation and animal diversity [J]. Nature, 2003, 424(6945):147-151.

[2] Smale T, Kadonaga T. The RNA polymerase II core promoter [J]. Ann Rev Biochem, 2003, 72(72):449-479.

[3] Qian J, Lin J, Luscombe N M, et al. Prediction of regulatory networks: genome-wide identification of transcription factor targets from gene expression data [J]. Bioinformatics, 2003, 19(15):1917-1926.

[4] Horsman S, Moorhouse M J, de Jager V C, et al. TF Target Mapper: a BLAST search tool for the identification of transcription factor target genes [J]. BMC Bioinformatics, 2006,7(1):1-7.

[5] Jolly E R, Chin C S, Herskowitz I, et al. Genome-wide identification of the regulatory targets of a transcription factor using biochemical characterization and computational genomic analysis [J]. BMC Bioinformatics, 2005,6(1):1-12.

[6] Christian D, Ralf S, Erich G, et al. MYB transcription factors inArabidopsis[J]. Cell, 2010,15(10):573-581.

[7] Riechmann J L, Ratcliffe O J. A genomic perspective on plant transcription factors [J]. Curr Opin Plant Biol, 2000,3(5):423-434.

[8] Stracke R, Ishihara H, Huep G, et al. Differential regulation of closely related R2R3-MYB transcription factors controls flavonol accumulation in different parts of theArabidopsisthalianaseedling [J]. Plant J, 2007,50(4):660-677.

[9] Fornale S, Sonbol F M, Capelledes M, et a1. ZmMYB31&ZmMYB42: two maize R2R3-MYB transcription factors having complementary roles in the liglin and phenylpropanoid metabolism regulation [J]. New Biotechnol, 2009, 25(S):5279-5280.

[10] Dugassa N F, Solveig M O, Behzad H. Nitrogen depletion and small R3-MYB transcription factors affecting anthocyanin accumulation inArabidopsisleaves [J]. Photochemistry, 2014,98(1):34-40.

[11] Zhao K, Bartley L E. Comparative genomic analysis of the R2R3 MYB secondary cell wall regulators ofArabidopsis, poplar, rice, maize, and switchgrass. BMC Plant Biol, 2014,14(1):135.

[12] Du H, Feng B R, Yang S S, et al. The R2R3-MYB transcription factor gene family in maize [J]. PLoS One, 2012,7(6):e37463.

[13] Song S S, Qi T C, Huang H, et al. The jasmonate-ZIM domain proteins interact with the R2R3-MYB transcription factors MYB21 and MYB24 to affect jasmonate-regulated stamen development inArabidopsis[J]. Plant cell, 2011,23(3):1000-1013.

[14] Velez-Bermudez I C, Salazar-Henao J E, Fornale S, et al. A MYB/ZML complex regulates wound-induced lignin genes in maize [J]. Plant Cell, 2015,27(11):3245-3259

[15] Blanchette M, Tompa M. Discovery of regulatory elements by a computational method for phylogenetic footprinting [J]. Genome Res, 2002,12(5):739-748.

[16] Huang W L, Tung C W, Huang H L, et al. ProLoc: prediction of protein subnuclear localization using SVM with automatic selection from physicochemical composition features [J]. Biosystems, 2007,90(2):573-581.

[17] Crooks G E, Hon G, Chandonia J M, et al. WebLogo: a sequence logo generator [J]. Genome Res, 2004,14(6):1188-1190.

[18] Du Z, Zhou X, Ling Y, et al. agriGO: a GO analysis toolkit for the agricultural community [J]. Nucleic Acid Res, 2010, 38:W64-70.

[19] Wang Q, Dooner H K. Remarkable variation in maize genome structure inferred from haplotype diversity at thebzlocus [J]. Proc Natl Acad Sci USA, 2006,103(47):17644-17649.

[20] Matys V, Fricke E, Geffers R, et al. TRANSFAC®: transcriptional regulation, from patterns to profiles [J]. Nucleic Acid Res, 2003,31(1):374-378.

[21] Bryne J C, Valen E, Tang M H E, et al. JASPAR, the open access database of transcription factor-binding profiles: new content and tools in the 2008 update [J]. Nucleic Acid Res, 2008, 36:D102-D106.