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

李琦 王怡超 劉暢 譚何新

（1. 海軍軍醫大學藥學院，上海 200433；2. 上海中醫藥大學岳陽中西醫結合醫院，上海 200437）

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

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

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

1 材料與方法

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

黃花蒿基因組數據庫、蛋白數據庫及注釋文件 下 載 于NCBI（National Center for Biotechnology Information），ID：PRJNA416223［15］， 利 用 NCBI BLAST 軟 件（ncbi-blast-2.9.0+-win64）構 建 本 地蛋白數據庫，從TAIR（The Arabidopsis Information Resource）數據庫（https://www.arabidopsis.org/）獲取擬南芥197條MYB轉錄因子序列［16］，并以此為查詢序列對數據庫進行本地 BLASTP 搜索，E-value值設為 0.0001；此外，在 Pfam 網站（https://Pfam.xfam.org/）下載 MYB 保守結構域的 HMM（Hidden Markov Model）模型文件 PF00249［17］，以此文件作為種子，應用 Hmmer 軟件，在本地蛋白數據庫中運行 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結構域及數量，根據結構特征篩選R2R3-MYB序列。使用 ProtParam在 線 工 具（https://web.expasy.org/protparam/）分析所得序列的分子量、等電點、GRAVY（grand average of hydropathicity）等理化信息［19］。

1.2 黃花蒿R2R3-MYB轉錄因子保守結構域分析及系統進化樹構建

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

1.3 黃花蒿R2R3-MYB轉錄因子功能預測

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

1.4 基于轉錄組的R2R3-MYB基因表達模式分析

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

2 結果

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

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

表1 黃花蒿R2R3-MYB轉錄因子篩選結果Table 1 Screening result of R2R3-MYB transcription factors in A. annua

利用ProtParam在線工具分析獲得的R2R3-MYB轉錄因子的分子量、等電點、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保守結構域分析

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

表2 黃花蒿R2R3-MYB轉錄因子MEME結構域搜索結果Table 2 Domain search result of R2R3-MYB transcription factors in Artemisia annua by MEME

圖1 黃花蒿R2R3-MYB轉錄因子保守結構域分析結果Fig.1 Conserved domain analysis result of R2R3-MYB transcription factors in A. annua

2.3 黃花蒿R2R3-MYB轉錄因子系統進化樹分析

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

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

2.4 黃花蒿R2R3-MYB轉錄因子功能注釋

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

表3 黃花蒿R2R3-MYB轉錄因子GO注釋結果統計Table 3 GO annotation count result of R2R3-MYB transcription factors in A. annua

2.5 基于轉錄組的R2R3-MYB基因表達模式分析

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

圖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轉錄因子基因，對其編碼的R2R3-MYB轉錄因子的理化性質進行了初步分析。對蛋白保守結構域的分析表明，黃花蒿R2R3-MYB轉錄因子均含有2個相鄰的MYB結構域，其結構符合“-W-（X19/20）-W-（X19）-W-…-F/I-（X18）-W-（X18）-W-”的特征，這表明黃花蒿與其它植物的R2R3-MYB轉錄因子在DNA結合域上具有較高的保守性［10］，而MEME軟件分析的結果除了能顯示DNA結合域的保守性之外，也能顯示出功能域部分的多樣性，這與R2R3-MYB轉錄因子功能廣泛相適應［29］。

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

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

4 結論

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