黄花蒿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转录因子进行了初步的筛选和分析，获得了可能参与青蒿素生物合成过程的功能基因，可为特定功能转录因子的筛选提供依据，其具体的生物学功能还需要验证。