李琦 王怡超 劉暢 譚何新

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

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

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

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

1 材料與方法

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

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

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

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

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

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

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

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

2 結(jié)果

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

構(gòu)建的黃花蒿本地蛋白數(shù)據(jù)庫(kù)，以擬南芥197條MYB轉(zhuǎn)錄因子的氨基酸序列為查詢序列，進(jìn)行本地BLASTP搜索，E-value值為0.0001的條件下共獲得470條序列；以HMM模型文件PF00249作為種子，使用HMMsearch搜索黃花蒿本地蛋白數(shù)據(jù)庫(kù)，E-value值為0.01的條件下共獲得597 條序列。合并所得序列，利用CD-Hit在線工具分析去除冗余序列，通過HMMscan在線分析確定序列中含有MYB結(jié)構(gòu)域的數(shù)量，篩選出符合R2R3-MYB轉(zhuǎn)錄因子特征的序列132條，編號(hào)為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)錄因子的分子量、等電點(diǎn)、GRAVY等信息。R2R3-MYB蛋白序列長(zhǎng)度為146-732個(gè)氨基酸殘基，平均309個(gè)；分子量17 270.05 Da-81 735.19 Da，平均35 067.82 Da；等電點(diǎn)4.56-9.82，平均6.96；GRAVY -1.068--0.359，平均-0.732。

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

Weblogo軟件對(duì)R2R3-MYB轉(zhuǎn)錄因子的保守結(jié)構(gòu)域進(jìn)行分析發(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)域形成相對(duì)保守的空間結(jié)構(gòu)。Swiss-Model軟件預(yù)測(cè)的三維結(jié)構(gòu)表明，R2和R3均形成了螺旋-螺旋-轉(zhuǎn)角-螺旋（helix-helix-turn-helix，HHTH）的結(jié)構(gòu)（圖1-C，D），與該區(qū)域結(jié)合DNA的功能相適應(yīng)［10，29］。此外，使用MEME在線軟件分析132條R2R3-MYB轉(zhuǎn)錄因子的保守基序，識(shí)別出的45個(gè)結(jié)構(gòu)域中E-value值最小的3個(gè)基序在幾乎所有序列中均出現(xiàn)，且都在R2R3結(jié)構(gòu)域中，表明黃花蒿R2R3-MYB轉(zhuǎn)錄因子的DNA結(jié)合域相對(duì)保守，而功能域則表現(xiàn)出了較大程度的多變性與復(fù)雜性（表2），這與R2R3-MYB轉(zhuǎn)錄因子功能廣泛相 適應(yīng)。

表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)進(jìn)化樹分析

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

圖2 黃花蒿與擬南芥R2R3-MYB轉(zhuǎn)錄因子系統(tǒng)進(jìn)化樹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蛋白共對(duì)應(yīng)了36個(gè)GO類別（表3）。其中在細(xì)胞組分（cellular component，C）分類中，有“細(xì)胞核”（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）分類中，有“細(xì)胞分化”（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)錄因子的功能相關(guān)。此外有些蛋白序列注釋有“響應(yīng)生長(zhǎng)素”（response to auxin，GO：0009733）、“氣孔運(yùn)動(dòng)調(diào)節(jié)”（regulation of stomatal movement，GO：0010090）等功能，這可能與該轉(zhuǎn)錄因子獨(dú)特的功能相關(guān)。

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

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

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

圖3 青蒿素特異合成途徑基因與部分黃花蒿R2R3-MYB基因表達(dá)模式聚類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)錄因子基因，對(duì)其編碼的R2R3-MYB轉(zhuǎn)錄因子的理化性質(zhì)進(jìn)行了初步分析。對(duì)蛋白保守結(jié)構(gòu)域的分析表明，黃花蒿R2R3-MYB轉(zhuǎn)錄因子均含有2個(gè)相鄰的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)錄因子功能廣泛相適應(yīng)［29］。

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

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

4 結(jié)論

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