玉米MYB轉(zhuǎn)錄因子靶基因的全基因組預測及驗證

張曉芳，于好強，羅 羲，張于福，李晚忱，付鳳玲

(四川農(nóng)業(yè)大學玉米研究所/農(nóng)業(yè)部西南玉米生物學與遺傳育種重點實驗室,四川 成都 611130)

轉(zhuǎn)錄因子對靶基因啟動子結(jié)合位點的識別與結(jié)合，是基因表達調(diào)控最關(guān)鍵的環(huán)節(jié)，也是功能基因組學分析的重要內(nèi)容[1,2]。以往通過酵母單雜交、染色質(zhì)免疫共沉淀(ChIP)、逆轉(zhuǎn)錄PCR(RT-PCR)、基因表達芯片等實驗手段篩選轉(zhuǎn)錄因子特異識別的啟動子，鑒別其靶基因，準確性不高，效率低。隨著生物信息學數(shù)據(jù)的積累，結(jié)合機器學習等智能計算技術(shù)的應用，根據(jù)轉(zhuǎn)錄因子結(jié)合位點以及轉(zhuǎn)錄因子與其靶基因表達模式的相關(guān)性，已開發(fā)出一些預測轉(zhuǎn)錄因子靶基因的生物信息學軟件[3-5]。

MYB是人類成髓細胞瘤(MYB)轉(zhuǎn)錄因子家族的植物直系同源蛋白，根據(jù)其包含MYB重復結(jié)構(gòu)的個數(shù)分為不同亞家族。在植物中，大部分MYB轉(zhuǎn)錄因子均包含R2和R3兩個MYB重復結(jié)構(gòu)，因此劃分為R2R3-MYB亞家族，調(diào)控植物眾多生長發(fā)育及逆境應答相關(guān)基因的表達[6-12]。玉米中已克隆的R2R3-MYB亞家族基因有C1、P1、MYB-IF25、MYB-IF35等，都與逆境應答有關(guān)[11-12]。但是，對于該亞家族更多MYB轉(zhuǎn)錄因子的靶基因及調(diào)控機制，卻鮮有報道[13]。R2R3-MYB亞家族轉(zhuǎn)錄因子識別啟動子的核心序列為TAACTG，其中第三個堿基A具有高度保守性，在MYB識別靶基因中起關(guān)鍵作用。但是，啟動子具有此核心序列的基因并不全是R2R3-MYB轉(zhuǎn)錄因子調(diào)控的靶基因。

為了鑒定更多R2R3-MYB轉(zhuǎn)錄因子調(diào)控的靶基因，進一步研究其在玉米逆境應答及生長發(fā)育過程中的功能，本研究根據(jù)MYB識別序列的核心序列及其側(cè)翼序列，結(jié)合運用HexDIFF算法和支持向量機(SVM)構(gòu)建分類模型，在玉米全基因組范圍內(nèi)對MYB轉(zhuǎn)錄因子的靶基因進行預測，根據(jù)功能注釋分析其可能的生物學功能。以電泳遷移率實驗(EMSA)在體外驗證MYB轉(zhuǎn)錄因子與預測序列的結(jié)合，并在玉米愈傷組織中瞬時表達檢測預測的MYB靶基因啟動子的活性。

1 材料和方法

1.1 模體框數(shù)據(jù)集構(gòu)建

綜合運用加州大學圣克魯茲分校開發(fā)的基因組瀏覽器UCSC Genome Browser(http://genome.ucsc.edu/)、真核啟動子數(shù)據(jù)庫EPD(http://www.epd.isb-sib.ch/)、轉(zhuǎn)錄因子數(shù)據(jù)庫TFD(http://www.ifti.org/)、轉(zhuǎn)錄調(diào)控區(qū)數(shù)據(jù)庫TRRD(http://wwwmgs.bionet.nsc.ru/mgs/gnw/trrd)、真核轉(zhuǎn)錄因子基因組結(jié)合位點和DNA結(jié)合譜數(shù)據(jù)庫TransFac(http://www.gene-regulation.com/)、轉(zhuǎn)錄因子與DNA結(jié)合位點模體數(shù)據(jù)庫Jaspar(http://jaspar.cgb.ki.se/cgi-bin/jaspar_db.pl)、轉(zhuǎn)錄起始位點數(shù)據(jù)庫dbTSS(http://dbtss.hgc.jp/)、植物順式作用調(diào)控元件及啟動子序列分析數(shù)據(jù)庫PlantCare(http://oberon.fvms.ugent.be:8080/PlantCARE/index.html)、植物順式調(diào)控DNA元件中的基序數(shù)據(jù)庫PLACE(http://www.dna.affrc.go.jp/htdocs/PLACE/)、植物啟動子數(shù)據(jù)庫PlantPromDB(http://mendel.cs.rhul.ac.uk/mendel.php?topic=plantprom)、美國國立生物技術(shù)信息中心文獻數(shù)據(jù)庫PubMed(www.ncbi.nlm.nih.gov/pubmed)等數(shù)據(jù)庫，搜索經(jīng)酵母單雜交、ChIP、RT-PCR、基因表達芯片等技術(shù)篩選鑒定的植物R2R3-MYB轉(zhuǎn)錄因子的靶基因。用PERL語言編寫腳本，分別從擬南芥資源網(wǎng)(www.arabidopsis.com)、禾本科資源網(wǎng)(www.gramene.org)和玉米基因組數(shù)據(jù)庫(www.maizesequence.org)下載轉(zhuǎn)錄起始位點上游2000 bp序列，用啟動子預測軟件Promoter 2.0(http://www.cbs.dtu.dk/services/Promoter/)預測其中的啟動子序列，分別構(gòu)建Arab+orysa_promoter和Zea_promoter數(shù)據(jù)庫，從中截取已驗證MYB轉(zhuǎn)錄因子結(jié)合位點核心序列TAACTG及其兩側(cè)各100 bp的序列，構(gòu)建陽性模體框數(shù)據(jù)集(Positive motif frame data)。從Arab+orysa_promoter數(shù)據(jù)庫中，檢索已驗證MYB轉(zhuǎn)錄因子以外但具有核心序列TAACTG的序列，截取其兩側(cè)共206 bp序列，構(gòu)建陰性模體框數(shù)據(jù)集(Negative motif frame data)；以TAACTG核心六聚體為中心，在兩側(cè)以A、T、C、G四種核苷酸為元素隨機函數(shù)生成100 bp的序列，構(gòu)建隨機模體框數(shù)據(jù)集(Random motif frame data)；再從Zea_promoter數(shù)據(jù)庫中，檢索核心序列TAACTG，截取其兩側(cè)共206 bp序列，構(gòu)建待測模體框數(shù)據(jù)集(Test motif frame data)。前3個數(shù)據(jù)集又合稱為訓練模體框數(shù)據(jù)集(Train motif frame data)。

1.2 玉米MYB轉(zhuǎn)錄因子識別啟動子預測

用PERL語言編寫腳本，在上述各模體框數(shù)據(jù)集序列中截取所有可能的6 bp六聚體，用HexDIFF算法分別計算其在陽性模體框數(shù)據(jù)集、陰性模體框數(shù)據(jù)集和隨機模體框數(shù)據(jù)集中出現(xiàn)的頻率R(h)[14]：

參照Huang等介紹的方法[16]，將陽性模體框數(shù)據(jù)集、陰性模體框數(shù)據(jù)集和隨機模體框數(shù)據(jù)集中所有可能六聚體的出現(xiàn)頻率，用智能計算的SVM，綜合使用線性和非線性方法分類。將訓練模體框數(shù)據(jù)集和待測模體框數(shù)據(jù)集歸一化處理，在一定的保證概率預測待測模體框數(shù)據(jù)集中R2R3-MYB轉(zhuǎn)錄因子識別的核心序列TAACTG及側(cè)翼序列，用WEBLOGO軟件(http://weblogo.berkeley.edu/logo.cgi)在線生成模體框堿基頻率圖[17]。

1.3 基因功能富集分析

用禾本科資源網(wǎng)(http://www.gramene.org)的GmeneMART在線軟件，選擇Ensembl Gene ID濾器，從Plant Genes的Zeamaysgenes(AGPv2)數(shù)據(jù)庫，導出所預測MYB轉(zhuǎn)錄因子靶基因的功能注釋。再將導出的靶基因功能注釋(GO)導入AgriGO在線軟件的Singular Enrichment Analysis(SEA)工具[18]，選擇玉米基因組ZeamaysAGPv2作背景，以0.05顯著水平、Fisher統(tǒng)計測試算法提交數(shù)據(jù)，得到靶基因的功能富集分析結(jié)果。

1.4 凝膠遷移率檢測

從玉米模式自交系B73的cDNA 樣品擴增R2R3-MYB轉(zhuǎn)錄因子基因ZmMYB-IF25的開放閱讀框序列，經(jīng)DNAworks在線軟件(https://omictools.com/dnaworks-tool)對密碼子偏好性優(yōu)化后，構(gòu)建原核表達載體pET-21a(+)-ZmMYB-IF25，轉(zhuǎn)化大腸桿菌BL21菌株，在0.4 mmol/L異丙基-β-D-硫代半乳糖苷(Isopropyl-β-D-thiogalactoside, IPTG)誘導下，37℃ 240 r/min 培養(yǎng)3 h，用12%十二烷基硫酸鈉聚丙烯酰胺凝膠電泳(Sodium dodecyl sulfate polyacrylamide gel electrophoresis, SDS-PAGE)進行表達檢測后，裝柱洗脫，分離純化MYB-IF25蛋白。

隨機選取30條預測模體框序列(表1)，生工生物工程(上海)股份有限公司合成，以已驗證的玉米靶基因Bz1和Bz2啟動子MYB核心序列TAACTG及兩側(cè)各106 bp序列為陽性對照，用地高辛標記后與純化的ZmMYB-IF25蛋白在5×EMSA結(jié)合緩沖液中混合，以2倍濃度的未標記雙鏈探針為冷競爭結(jié)合對照，室溫保持20 min上樣，進行連續(xù)非變形聚丙烯酰胺凝膠電泳，轉(zhuǎn)膜和紫外交聯(lián)后化學發(fā)光并成像，根據(jù)顯色條帶的遷移率鑒別MYB-IF25蛋白與預測模體框序列的結(jié)合。

1.5 瞬時表達驗證

從體外驗證的預測模體框序列中隨機選取4條，擴增其全長啟動子序列，構(gòu)建成啟動β-葡萄糖苷酸酶(β-glucuronidase)基因GUS的瞬時表達載體pBI221-pCandidate-GUS，混合1/3熒光素酶(Luciferase)基因LUC表達載體pUbi-LUC為內(nèi)參，用DuPont PDS1000/He型基因槍(Bio-Rad, USA)，在1100 psi系統(tǒng)壓力和25 mm汞柱真空度下，轉(zhuǎn)化玉米胚性愈傷組織。每個處理重復3次，每次重復轉(zhuǎn)化15塊大小、形狀和質(zhì)地基本一致的愈傷組織。轉(zhuǎn)化后在8%甘露醇高滲培養(yǎng)基上27℃黑暗培養(yǎng)24 h，GUS組織化學染色后在解剖鏡下觀察拍照藍色斑點。剩余愈傷組織用熒光素酶檢測試劑盒(Promega，上海)顯色后，用Fluoroskan Ascent FL型熒光/化學發(fā)光微孔檢測儀(天齊生物科技有限公司，上海)，檢測高滲脅迫0(對照)和24 h的GUS與LUC熒光強度比值GUS/LUC，對預測靶序列響應高滲脅迫的體內(nèi)啟動活性進行相對定量。

表1 隨機選取的預測模體框序列

2 結(jié)果與分析

2.1 玉米MYB靶基因及啟動子結(jié)合位點

當SVM參數(shù)C=2，γ=0.03125時，以98.38%的保證概率，從Zea_promoter數(shù)據(jù)庫136770條序列中，預測到435個MYB轉(zhuǎn)錄因子結(jié)合位點，對應下游424個靶基因。其中，Bz1(GRMZM2G165390)和Bz2(GRMZM2G016241)是已驗證的MYB轉(zhuǎn)錄因子靶基因[19]。與玉米基因組序列比對結(jié)果，絕大部分預測靶基因隨機分布于全部10個玉米染色體上，沒有染色體偏好，但有的區(qū)段有集中分部情況，沒有分布于線粒體和葉綠體基因組，極少部分預測靶基因不能與基因組序列匹配(圖1)。

預測靶基因啟動子MYB識別位點及其側(cè)翼序列的堿基頻率與已知MYB靶基因相似，具有高度保守性，但距離核心序列較遠的堿基存在一定的變異(圖2)。這既與物種差異有關(guān)，也可能是同一物種基因間的差異。

圖1 預測MYB靶基因在玉米染色體上的分布Fig.1 Distribution of predicted MYB targeted genes on the maize chromosomes

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

用GrameneMART軟件從Plant Genes的Zeamaysgenes(AGPv2)數(shù)據(jù)庫中共搜索到229個預測MYB靶基因的功能注釋共計1471條。其中，45個(19.7%)基因涉及刺激響應(GO：0050896)，33個(14.4%)基因介導逆境應答(GO：0006950)，10個(4.3%)基因與再生過程(GO：0000003)相關(guān)，10個(4.3%)基因參與發(fā)育過程(GO：0032502)。以細胞組分類分析發(fā)現(xiàn)，93個(40.6%)基因編碼蛋白定位于細胞內(nèi)(GO：0005622)，12個(5.2%)基因編碼蛋白定位于細胞外(GO：0005576)。

2.3 MYB轉(zhuǎn)錄因子與預測靶基因啟動子序列的體外結(jié)合

EMSA結(jié)果表明，隨機選取的30條預測模體框序列中，除基因ID為GRMZM2G086773、GRMZM2G137596和GRMZM2G087719的3個基因的預測模體框序列以外(表1)，其余27條(90%)預測模體框序列均可在體外與MYB-IF25蛋白結(jié)合，電泳遷移條帶滯后(圖3)，表明本研究所用的轉(zhuǎn)錄因子全基因組預測方法具有較高的可靠性。

圖2 MYB靶基因核心及側(cè)翼序列堿基頻率Fig.2 Base frequency during core and flanking sequences of MYB targeted genes

FP: 無MYB-IF25蛋白結(jié)合的自由探針；B：MYB-IF25蛋白與地高辛標記預測模體框序列結(jié)合；C：MYB-IF25蛋白與未標記預測模體框序列結(jié)合(冷競爭結(jié)合)；1～30：隨機選取的30條預測模體框序列；bz1和bz2：Bz1和Bz2基因核心序列TAACTG及兩側(cè)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蛋白與預測模體框序列結(jié)合的EMSA檢測Fig.3 EMSA detection for combination between the ZmMYB-IF25 protein and the predicted motif frame sequences

2.4 預測靶基因啟動子序列的體內(nèi)表達活性

從體外驗證的27條預測模體框序列中隨機選取4條涉及逆境刺激響應相關(guān)靶基因(GRMZM2G044829、GRMZM2G115698、GRMZM2G310161和GRMZM2G430675)，擴增其全長啟動子序列，構(gòu)建成啟動GUS基因的瞬時表達載體pBI221-pCandidate-GUS，基因F槍法轉(zhuǎn)化玉米愈傷組織，8%甘露醇高滲培養(yǎng)基27℃黑暗培養(yǎng)24 h后，GUS組織化學染色結(jié)果顯示，4個預測靶基因的啟動子均能在玉米愈傷組織內(nèi)啟動GUS基因的表達，產(chǎn)生顯色反應(圖4)。這4個啟動子8%甘露醇滲透脅迫下啟動GUS表達產(chǎn)生的GUS熒光強度與LUC熒光強度的比值，均比脅迫前顯著增高(圖5)，說明這4個靶基因的啟動子在玉米細胞內(nèi)具有啟動活性，且受高滲脅迫誘導。

3 結(jié)論與討論

本研究將HexDIFF算法與SVM結(jié)合，用已知MYB結(jié)合位點序列建模，在玉米全基因組范圍內(nèi)預測到424個MYB靶基因和435個MYB結(jié)合位點(圖1、2)，涉及逆境刺激響應、再生和發(fā)育等眾多生長發(fā)育過程。EMSA實驗表明，預測的結(jié)合位點與玉米MYB可相互結(jié)合(圖3)。GUS瞬時表達證實預測得到的MYB靶基因啟動子具有啟動活性(圖4、5)，而且預測的逆境相關(guān)靶基因啟動子啟動報告基因GUS，在滲透脅迫下的表達量顯著提高，更說明預測的MYB靶基因有一部分的功能確與植物抗逆反應相關(guān)。以上結(jié)果證明本研究預測方法可靠性高，為預測轉(zhuǎn)錄因子靶基因提供借鑒。隨著第三代測序技術(shù)(尤其是轉(zhuǎn)錄組學)的快速發(fā)展，人們必將對玉米MYB轉(zhuǎn)錄因子進行更為深入的研究，屆時，玉米MYB轉(zhuǎn)錄因子參與的生長、發(fā)育等調(diào)控機制將會被更好地解析。

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體外驗證實驗發(fā)現(xiàn)，預測結(jié)果可能存在10%假陽性，可能的原因是各數(shù)據(jù)庫能夠檢索到的MYB靶基因及其結(jié)合位點有限，構(gòu)建的陽性模體框數(shù)據(jù)集代表性可能不夠全面。隨著轉(zhuǎn)錄因子靶基因數(shù)據(jù)量的積累，本研究方法可能會得到更好的驗證。本研究的分析還表明，核心序列TAACTG是MYB結(jié)合位點的必要條件，在我們檢索的陽性數(shù)據(jù)和預測的全部玉米MYB結(jié)合位點中100%保守，但并不是充分條件，MYB的的結(jié)合還決定于核心序列的側(cè)翼序列(圖2)。用TRANSFAC(http://gene-regulation.com/pub/databases.html)和JASPAR(http://jaspar.genereg.net/)等軟件搜索可發(fā)現(xiàn)， TAACTG序列還存在于DREB(Dehydration responsive element binding)、NAC(NAM/ATAF1/2/CUC2)等轉(zhuǎn)錄因子的啟動子序列中[20-21]。這也說明，MYB對其靶基因的轉(zhuǎn)錄調(diào)控可能還與DREB、NAC等轉(zhuǎn)錄因子存在協(xié)同關(guān)系。

參考文獻：

[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.