黃花蒿R2R3-MYB轉(zhuǎn)錄因子全基因組鑒定及生物信息學分析

李琦 王怡超 劉暢 譚何新

（1. 海軍軍醫(yī)大學藥學院，上海 200433；2. 上海中醫(yī)藥大學岳陽中西醫(yī)結(jié)合醫(yī)院，上海 200437）

以青蒿素為基礎的聯(lián)合用藥（artemisinin-based combination therapy，ACT）是世衛(wèi)組織推薦的瘧疾一線用藥［1］。青蒿素及其衍生物還被用于癌癥、紅斑狼瘡、風濕等疾病的治療中［2］，應用的擴大導致需求量進一步增大。因此，提升青蒿素產(chǎn)量是當前國際研究熱點。目前，青蒿素的唯一天然來源是藥用植物黃花蒿（Artemisia annua L.），但是其含量相對比較低，只占干重的0.1%-0.8%［3］。近年來，合成生物學在青蒿素生物半合成上取得了比較大的發(fā)展，已可較高效的合成青蒿酸（artemisinic acid）［4］，然而酵母無法像黃花蒿腺毛一樣提供青蒿素合成所需的特殊油性氧化環(huán)境，難以實現(xiàn)青蒿素的體外生物全合成。目前青蒿素市場供應仍然依賴于從黃花蒿中提取，培育高品質(zhì)黃花蒿株系意義重大［5-6］。

植物次生代謝過程中，轉(zhuǎn)錄因子可調(diào)控代謝途徑中一系列基因的表達，對轉(zhuǎn)錄因子的干預是一種有效的調(diào)控植物次生代謝產(chǎn)物合成的方法［7］，黃花蒿中已有多個家族的轉(zhuǎn)錄因子被證實參與青蒿素的生物合成途徑［8］。MYB（v-myb avian myeloblastosis viral oncogene homolog）轉(zhuǎn)錄因子家族具有保守的MYB結(jié)構(gòu)域（MYB domain），數(shù)量龐大且功能多樣，廣泛存在于所有真核生物中［9-10］。植物中MYB轉(zhuǎn)錄因子的保守結(jié)構(gòu)域通常由1-4個相鄰的重復序列組成，根據(jù)重復序列的數(shù)量和位置不同，MYB轉(zhuǎn)錄因子可分為4個亞家族，即R3-MYB、R2R3-MYB、R1R2R3-MYB和atypical-MYB，其中R2R3-MYB是最常見的類型［10］。MYB轉(zhuǎn)錄因子廣泛參與植物重要的生物過程，包括應激和防御反應、初生代謝與次生代謝的調(diào)控、細胞分化、種子和花的發(fā)育等［11］。黃花蒿中已有研究表明R2R3-MYB轉(zhuǎn)錄因子在青蒿素生物合成中發(fā)揮重要的調(diào)控作用，如AaMYB1可以調(diào)控青蒿素合成途徑基因的表達水平并增加分泌性腺毛數(shù)量及密度，從而影響青蒿素的合成與積累［12］；AaMIXTA1可以調(diào)控分泌性腺毛的數(shù)量和角質(zhì)層的生物合成，過表達AaMIXTA1可以在不影響分泌性腺毛結(jié)構(gòu)的情況下增加青蒿素的生成量［13］；AaTAR2能夠影響分泌性腺毛的形成和發(fā)育，調(diào)控青蒿素和黃酮類物質(zhì)的生物合成過程［14］。

相對于植物中所蘊含的龐大的數(shù)量，已報道的黃花蒿MYB轉(zhuǎn)錄因子個數(shù)僅僅是冰山一角，目前的研究仍不足，許多有潛力的轉(zhuǎn)錄因子等待挖掘，尤其是數(shù)量最多、功能最多樣的R2R3-MYB轉(zhuǎn)錄因子。Shen等［15］于2018年發(fā)布了黃花蒿基因組序列，為黃花蒿的研究提供了重要的參考數(shù)據(jù)。本研究基于已報道黃花蒿基因組、轉(zhuǎn)錄組數(shù)據(jù)，利用生物信息學分析的方法，在基因組層面上挖掘花蒿R2R3-MYB轉(zhuǎn)錄因子基因，為篩選調(diào)控青蒿素生物合成等重要生物過程的功能基因提供參考。

1 材料與方法

1.1 黃花蒿R2R3-MYB蛋白序列的獲取

黃花蒿基因組數(shù)據(jù)庫、蛋白數(shù)據(jù)庫及注釋文件 下 載 于NCBI（National Center for Biotechnology Information），ID：PRJNA416223［15］， 利 用 NCBI BLAST 軟 件（ncbi-blast-2.9.0+-win64）構(gòu) 建 本 地蛋白數(shù)據(jù)庫，從TAIR（The Arabidopsis Information Resource）數(shù)據(jù)庫（https://www.arabidopsis.org/）獲取擬南芥197條MYB轉(zhuǎn)錄因子序列［16］，并以此為查詢序列對數(shù)據(jù)庫進行本地 BLASTP 搜索，E-value值設為 0.0001；此外，在 Pfam 網(wǎng)站（https://Pfam.xfam.org/）下載 MYB 保守結(jié)構(gòu)域的 HMM（Hidden Markov Model）模型文件 PF00249［17］，以此文件作為種子，應用 Hmmer 軟件，在本地蛋白數(shù)據(jù)庫中運行 HMMsearch，E-value值設為0.01。合并所得序列，使用CD-HIT在線分析（https://www.bioinformiscs.org/CD-HIT/）去除冗余序列［18］，對剩余序列進行 HMMscan 在 線 分 析（https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan），確定MYB結(jié)構(gòu)域及數(shù)量，根據(jù)結(jié)構(gòu)特征篩選R2R3-MYB序列。使用 ProtParam在 線 工 具（https://web.expasy.org/protparam/）分析所得序列的分子量、等電點、GRAVY（grand average of hydropathicity）等理化信息［19］。

1.2 黃花蒿R2R3-MYB轉(zhuǎn)錄因子保守結(jié)構(gòu)域分析及系統(tǒng)進化樹構(gòu)建

在Clustal X2軟件對獲得的R2R3-MYB轉(zhuǎn)錄因子序列進行對齊，截取保守結(jié)構(gòu)域，通過 Weblogo 網(wǎng)站（https://weblogo.berkeley.edu/）進行保守氨基酸基序分析［20］，并用Swiss-Model在線軟件（https://swissmodel.expasy.org/）預測保守結(jié)構(gòu)域的三維結(jié)構(gòu)［21］；此外利用MEME5.1.1在線軟件（https://meme-suite.org）分析R2R3-MYB轉(zhuǎn)錄因子的保守基序［22］，根據(jù)文獻報道，將最大結(jié)構(gòu)域數(shù)量設置為45，其余使用默認參數(shù)［23］。以擬南芥126條 R2R3-MYB轉(zhuǎn)錄因子序列為參考，利用Mega6.0軟件進行序列對齊，選擇鄰近法（Neighbor-Joining）構(gòu)建黃花蒿與擬南芥R2R3-MYB 轉(zhuǎn)錄因子的系統(tǒng)進化樹［24］。

1.3 黃花蒿R2R3-MYB轉(zhuǎn)錄因子功能預測

通過Blast2GO v5.2.5軟件（https://www.blast2go.com/）對黃花蒿MYB轉(zhuǎn)錄因子的功能進行預測［25］，使用軟件中的Blast工具，以黃花蒿R2R3-MYB轉(zhuǎn)錄因子序列在NCBI非冗余數(shù)據(jù)庫（Nr Database）搜索相似序列，E-value值設為10-3，比對完成后使用默認參數(shù)在Blast2GO上執(zhí)行映射和注釋。

1.4 基于轉(zhuǎn)錄組的R2R3-MYB基因表達模式分析

黃花蒿轉(zhuǎn)錄組數(shù)據(jù)下載于NCBI（PRJNA416223［15］，PRJNA39657［26］），使 用NCBI SRA Toolkit軟件（version 2.10.0）解壓縮。黃花蒿mRNA序列由TBtools軟件從基因組中提取獲得［27］，并用Salmon軟件構(gòu)建索引，對轉(zhuǎn)錄組數(shù)據(jù)比對計數(shù)后計算TPM值（Transcripts PerKilobase Million）［28］，計算結(jié)果導入MEV4.9.0軟件制作熱圖并進行聚類分析。

2 結(jié)果

2.1 黃花蒿R2R3-MYB蛋白序列的獲取及理化性質(zhì)分析

構(gòu)建的黃花蒿本地蛋白數(shù)據(jù)庫，以擬南芥197條MYB轉(zhuǎn)錄因子的氨基酸序列為查詢序列，進行本地BLASTP搜索，E-value值為0.0001的條件下共獲得470條序列；以HMM模型文件PF00249作為種子，使用HMMsearch搜索黃花蒿本地蛋白數(shù)據(jù)庫，E-value值為0.01的條件下共獲得597 條序列。合并所得序列，利用CD-Hit在線工具分析去除冗余序列，通過HMMscan在線分析確定序列中含有MYB結(jié)構(gòu)域的數(shù)量，篩選出符合R2R3-MYB轉(zhuǎn)錄因子特征的序列132條，編號為AaMYB2R1到AaMYB2R132（表1）。

表1 黃花蒿R2R3-MYB轉(zhuǎn)錄因子篩選結(jié)果Table 1 Screening result of R2R3-MYB transcription factors in A. annua

利用ProtParam在線工具分析獲得的R2R3-MYB轉(zhuǎn)錄因子的分子量、等電點、GRAVY等信息。R2R3-MYB蛋白序列長度為146-732個氨基酸殘基，平均309個；分子量17 270.05 Da-81 735.19 Da，平均35 067.82 Da；等電點4.56-9.82，平均6.96；GRAVY -1.068--0.359，平均-0.732。

2.2 黃花蒿R2R3-MYB保守結(jié)構(gòu)域分析

Weblogo軟件對R2R3-MYB轉(zhuǎn)錄因子的保守結(jié)構(gòu)域進行分析發(fā)現(xiàn)，R2結(jié)構(gòu)域在4位、25位、45位，R3結(jié)構(gòu)域的78位、97位為保守的色氨酸（W）殘基，符合 “-W-（X19/20）-W-（X19）-W-…-F/I-（X18）-W-（X18）-W-”的結(jié)構(gòu)特征［10］，此外還有一些氨基酸具有高度的保守性，如谷氨酸（E）、天冬氨酸（D）、亮氨酸（L）、精氨酸（R）、賴氨酸（K）、半胱氨酸（C）、甘氨酸（G）、天冬酰胺（N）等（圖1-A，B），這使R2R3結(jié)構(gòu)域形成相對保守的空間結(jié)構(gòu)。Swiss-Model軟件預測的三維結(jié)構(gòu)表明，R2和R3均形成了螺旋-螺旋-轉(zhuǎn)角-螺旋（helix-helix-turn-helix，HHTH）的結(jié)構(gòu)（圖1-C，D），與該區(qū)域結(jié)合DNA的功能相適應［10，29］。此外，使用MEME在線軟件分析132條R2R3-MYB轉(zhuǎn)錄因子的保守基序，識別出的45個結(jié)構(gòu)域中E-value值最小的3個基序在幾乎所有序列中均出現(xiàn)，且都在R2R3結(jié)構(gòu)域中，表明黃花蒿R2R3-MYB轉(zhuǎn)錄因子的DNA結(jié)合域相對保守，而功能域則表現(xiàn)出了較大程度的多變性與復雜性（表2），這與R2R3-MYB轉(zhuǎn)錄因子功能廣泛相 適應。

表2 黃花蒿R2R3-MYB轉(zhuǎn)錄因子MEME結(jié)構(gòu)域搜索結(jié)果Table 2 Domain search result of R2R3-MYB transcription factors in Artemisia annua by MEME

圖1 黃花蒿R2R3-MYB轉(zhuǎn)錄因子保守結(jié)構(gòu)域分析結(jié)果Fig.1 Conserved domain analysis result of R2R3-MYB transcription factors in A. annua

2.3 黃花蒿R2R3-MYB轉(zhuǎn)錄因子系統(tǒng)進化樹分析

利用Mega6.0軟件構(gòu)建黃花蒿與擬南芥中R2R3-MYB轉(zhuǎn)錄因子的系統(tǒng)進化樹，將黃花蒿的132條序列與擬南芥的126條序列對齊后，分析構(gòu)建系統(tǒng)進化樹的最佳模型，根據(jù)軟件計算所得的BIC值（bayesian information criterion），選擇分值最小的JTT模型（Jones-Taylor-Thornton model），采用鄰近法（neighbor-joining）構(gòu)建系統(tǒng)進化樹，并根據(jù)已有的擬南芥文獻報道，將黃花蒿R2R3-MYB分入23個亞組［10，29］，未分入亞組的按G1-G18編號分組（圖2）。通過分析進化樹發(fā)現(xiàn)，擬南芥的S3、S12、S15亞組中沒有黃花蒿序列，一些組中也沒有擬南芥序列，如G1、G2、G5、G6、G8、G16，這也反映了植物R2R3-MYB轉(zhuǎn)錄因子在進化過程中的復雜性和多樣性。

圖2 黃花蒿與擬南芥R2R3-MYB轉(zhuǎn)錄因子系統(tǒng)進化樹Fig.2 Phylogenetic tree of R2R3-MYB transcription factors in A. annua and Arabidopsis thaliana

2.4 黃花蒿R2R3-MYB轉(zhuǎn)錄因子功能注釋

通過整理Blast2-GO軟件分析結(jié)果，本研究中篩選得到的132條R2R3-MYB蛋白共對應了36個GO類別（表3）。其中在細胞組分（cellular component，C）分類中，有“細胞核”（nucleus，GO：0003677）注釋的蛋白數(shù)量為128，占97%。在分子功能（molecular function，F(xiàn)）分類中，有“特定序列DNA結(jié)合”（sequence-specific DNA binding，GO：0043565）注釋的數(shù)量為95，占72.0%；有“轉(zhuǎn)錄調(diào)控區(qū)DNA結(jié)合”（transcription regulatory region DNA binding，GO：0044212）注釋的數(shù)量為84，占63.6%；有“DNA結(jié)合轉(zhuǎn)錄因子活性”（DNA-binding transcription factor activity，GO：0003700）注釋的數(shù)量為83，占62.9%；有“結(jié)合DNA”（DNA binding，GO：0003677）注釋的數(shù)量為36，占27.3%。在生物過程（biological process，P）分類中，有“細胞分化”（cell differentiation，GO：0030154）注釋的數(shù)量為80，占60.6%；有“DNA為模板的轉(zhuǎn)錄調(diào)控”（regulation of transcription，DNA-templated，GO：0006355）注釋的數(shù)量為67，占50.8%。這些注釋均與R2R3-MYB蛋白作為轉(zhuǎn)錄因子的功能相關。此外有些蛋白序列注釋有“響應生長素”（response to auxin，GO：0009733）、“氣孔運動調(diào)節(jié)”（regulation of stomatal movement，GO：0010090）等功能，這可能與該轉(zhuǎn)錄因子獨特的功能相關。

表3 黃花蒿R2R3-MYB轉(zhuǎn)錄因子GO注釋結(jié)果統(tǒng)計Table 3 GO annotation count result of R2R3-MYB transcription factors in A. annua

2.5 基于轉(zhuǎn)錄組的R2R3-MYB基因表達模式分析

本研究利用Salmon軟件計算獲得了R2R3-MYB基因和青蒿素特異合成途徑酶基因（ADS、CYP71AV1、DBR2、ALDH1）在不同組織表達的TPM值，其中PRJNA416223轉(zhuǎn)錄組數(shù)據(jù)可分析基因在根、莖、幼葉、老葉、表皮、花蕾、花、種子、腺毛9個組織中的表達水平（圖3-A），PRJNA39657轉(zhuǎn)錄組數(shù)據(jù)可分析基因在分生組織、幼葉腺毛、成葉腺毛、花蕾腺毛和子葉中的表達水平（圖3-B）。根據(jù)計算得到的TPM值使用MEV4.9.0軟件繪制熱圖并進行聚類分析，在兩個轉(zhuǎn)錄組的分析結(jié)果中，MYB2R60（mRNA ID：AA279750.t1）基 因 與青蒿素特異合成途徑酶基因的表達模式均相近，并被聚類在一起，這與基因組文章中預測的結(jié)果一致［15］；同 時MYB2R27（mRNA ID：AA156700.t1）和MYB2R96（即AaMYB1）［12］基因在不同組織腺毛中與CYP71AV1基因有相近的表達模式并被聚類在一起，可能參與到CYP71AV1基因表達的調(diào)控。

圖3 青蒿素特異合成途徑基因與部分黃花蒿R2R3-MYB基因表達模式聚類Fig.3 Expression pattern clustering of artemisinin specific synthesis pathway genes and parts of R2R3-MYB genes in A. annua

3 討論

本研究在基因組水平上篩選獲得了132條黃花蒿R2R3-MYB轉(zhuǎn)錄因子基因，對其編碼的R2R3-MYB轉(zhuǎn)錄因子的理化性質(zhì)進行了初步分析。對蛋白保守結(jié)構(gòu)域的分析表明，黃花蒿R2R3-MYB轉(zhuǎn)錄因子均含有2個相鄰的MYB結(jié)構(gòu)域，其結(jié)構(gòu)符合“-W-（X19/20）-W-（X19）-W-…-F/I-（X18）-W-（X18）-W-”的特征，這表明黃花蒿與其它植物的R2R3-MYB轉(zhuǎn)錄因子在DNA結(jié)合域上具有較高的保守性［10］，而MEME軟件分析的結(jié)果除了能顯示DNA結(jié)合域的保守性之外，也能顯示出功能域部分的多樣性，這與R2R3-MYB轉(zhuǎn)錄因子功能廣泛相適應［29］。

對系統(tǒng)進化樹進行分析表明，部分黃花蒿R2R3-MYB轉(zhuǎn)錄因子序列與擬南芥序列同源性較高，或可分入同一亞組，提示可能與對應的擬南芥轉(zhuǎn)錄因子具有相似的生物功能，如AaMYB2R30（即AaMIXTA1）［13］與擬南芥S9亞組轉(zhuǎn)錄因子高度同源、AaMYB2R96（即AaMYB1）［12］與擬南芥S13亞組轉(zhuǎn)錄因子高度同源；同時，部分黃花蒿R2R3-MYB轉(zhuǎn)錄因子所在的組中不含有擬南芥序列，這些黃花蒿R2R3-MYB轉(zhuǎn)錄因子則可能具有相對于擬南芥而言獨特的生物學功能。

在對黃花蒿R2R3-MYB蛋白功能的預測中，絕大多數(shù)蛋白均被預測定位于細胞核中，且功能上均與結(jié)合DNA（DNA binding）相關，這是R2R3-MYB轉(zhuǎn)錄因子作為轉(zhuǎn)錄因子的共性特征。此外，部分序列擁有獨特的功能注釋，如AaMYB2R24（即TAR2）［14］有“響應赤霉素”（response to gibberellin，GO：0009739）、“響應水楊酸”（response to salicylic acid，GO：0009751）和“響應茉莉酸”（response to jasmonic acid，GO：0009753）的注釋，提示該轉(zhuǎn)錄因子可能參與到黃花蒿對激素響應的途徑中。同時本研究借助公共轉(zhuǎn)錄組數(shù)據(jù)，通過軟件計算基因的表達量，預測基因在不同組織中的表達模式，通過與青蒿素特異合成途徑中關鍵基因的表達模式比對，篩選到可能與青蒿素合成相關的基因，如AaMYB1［12］即本研究中的Aa2RMYB96，計算獲得的數(shù)據(jù)及比對方法也可以用于篩選其他生物途徑中的相關基因。

4 結(jié)論

本研究使用生物信息學的方法，對黃花蒿R2R3-MYB轉(zhuǎn)錄因子進行了初步的篩選和分析，獲得了可能參與青蒿素生物合成過程的功能基因，可為特定功能轉(zhuǎn)錄因子的篩選提供依據(jù)，其具體的生物學功能還需要驗證。