Genome-wide identification and expression pattern analysis of DoMYB gene family in Dalbergia odorifera

3年生降香黄檀植株样本采自海南省海南大学儋州校区试验基地（19°30'40" N，109°29'25" E），取样时间为2021年11月；7年生降香黄檀植株样本采自海南省定安县人工培育苗圃（19°38'56" N，110°14'29" E），取样时间为2020年7月。分别采集3年生与7年生降香黄檀植株的韧皮部、边材、茎干过渡区及叶片组织，用于转录组测序（RNA-seq）分析，其中，茎干过渡区定义为植株茎干心材与边材之间的过渡生长区域。本氏烟草（Nicotiana benthamiana）作为亚细胞定位试验材料，常规培养条件设置为：温度24℃、光照周期16 h光照/8 h黑暗。

本研究所用降香黄檀基因组数据为本课题组前期组装注释的未公开高质量基因组，基因组总大小约 1.2 Gb，Contig N50 为 58.21 Mb，共注释得到56 679 个蛋白编码基因，基因组组装与注释完整性良好。

从 Pfam 数据库（https://pfam.xfam.org/）获取 MYB 保守结构域对应的隐马尔可夫模型 HMM 文件 PF00249 ^[17]，借助 TBtools 软件中 Simple HMM Search 工具对降香黄檀全基因组蛋白序列进行比对检索，筛选标准设定 E-value＜1×10⁻⁵。同时，从拟南芥专用数据库 TAIR（http://www.arabidopsis.org/）下载拟南芥全部 MYB 家族蛋白序列作为查询集 ^[18]，采用 BLASTP 程序比对降香黄檀蛋白库，阈值同样设为 E-value＜1×10⁻⁵。

合并 HMM 检索与 BLASTP 比对两套筛选结果，剔除重复冗余序列；再分别利用 Pfam 在线工具与 NCBI 保守结构域批量检索平台 Batch CD-search（https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi?SEQUENCE）校验候选蛋白的 MYB 结构域完整度，删除结构域残缺的序列，最终获得可靠的 DoMYB 家族成员 ^[19]。

采用 Expasy ProtParam 工具（https://web.expasy.org/protparam/）批量预测 DoMYB 蛋白的分子量、理论等电点、亲疏水性、不稳定系数等基础理化特征 ^[20]；通过 WoLF PSORT 在线平台（https://wolfpsort.hgc.jp/）完成所有 MYB 蛋白的亚细胞定位生物信息学预测，物种分类选择植物（Plant）模式开展分析 ^[21]。

根据基因组注释文件，利用TBtools软件将降香黄檀DoMYB基因家族成员在染色体上的物理位置进行可视化展示，并对基因结构进行分析。使用NCBI Batch CD-search进行保守结构域的预测；运用MEME Suite（https://meme-suite.org/meme/）在线分析工具对保守基序进行预测，基序识别的最大数量被设定为20。使用TBtools进行保守结构域、保守基序及基因结构的可视化分析 ^[22]。

从TAIR数据库、Uniprot数据库（https://www.uniprot.org/）和Phytozome数据库（https://phytozome-next.jgi.doe.gov/）下载拟南芥、大豆（Glycine max）、水稻、毛果杨的家族的数据。利用ClustalW算法对拟南芥与降香黄檀的序列进行比对，使用MEGA12（v12.0）软件基于邻接法（Neighbor-Joining，NJ）构建系统发育树，bootstrap重复次数为1 000次，其余参数为默认值。降香黄檀种内系统进化树划分为23个亚组（A1～A23），与拟南芥MYB基因构建的系统发育树划分为25个亚组（C1～C25）。使用iTOL（https://itol.embl.de/）进行美化和可视化。在生成的系统发育树中，拟南芥的基因用方框标记，降香黄檀的基因用星形标记。此外，不同的MYB亚家族成员用颜色清晰区分：1R-MYB蛋白以棕色标示，2R-MYB蛋白以黑色标示，3R-MYB蛋白以红色标示，拟南芥中非典型的4R-MYB蛋白则以蓝色标示。

使用TBtools中的MCScanX插件对降香黄檀种内的MYB基因进行共线性分析，同时计算非同义替换率（Ka）、同义替换率（Ks）以及二者的比值（Ka/Ks），其中Ka/Ks>1表明正选择，Ka/Ks=1表明中性进化，Ka/Ks<1表明纯化选择；使用TBtools中的MCScanX插件对降香黄檀（Do）、拟南芥（At）、水稻（Os）、大豆（Gm）、毛果杨（Ptr）的基因组内的MYB基因进行共线性分析，绘制物种间共线性图谱并可视化展示共线性基因对。

本研究采用TBtools程序，对所有转录起始位点上游2 000 bp范围内的启动子片段进行批量提取，再将得到的序列导入PlantCARE（https://bioinformatics.psb.ugent.be/webtools/plantcare/html）在线网站进行顺式调控元件的识别与注释，对预测结果进行整理筛选与精简归类，使用Excel进行绘制。

利用降香黄檀3年生树木的韧皮部（B01～B05）、边材（S01～S05）、过渡区（T01～T05）及叶片（L01～L05）四个组织的转录组表达数据（PRJCA056155）的TPM值，其中每个组织的样本取自5株已形成心材的独立植株，使用Z-score标准化数据进行热图绘制。采用RNA提取试剂盒（CW0581S），从降香黄檀7年生植株的韧皮部、边材、过渡区和叶片组织中提取RNA。利用1 μg RNA，按照PrimeScript TM RT reagent Kit with gDNA Eraser试剂盒配置反转录体系，反转录为cDNA；以TB Green Premix Ex Taq Ⅱ试剂盒配置qPCR体系，设计DoMYB48基因特异性引物，以Actin为内参基因进行RT-qPCR反应 ^[23]。实验设置3次生物学重复，每次重复包含3次技术重复。其中RT-qPCR反应体系如下：引物为0.4 μL、TB Green Premix Ex TaqII为10 μL、cDNA为2 μL、无菌水为7.2 μL；反应程序使用两步法PCR扩增程序：95℃预变性30 s，95℃条件下5 s，60℃条件下34 s，循环数为40。使用2^-∆∆Ct法计算相对表达量，使用Graphpad Prism制作柱状图，采用单因素方差（One-way ANOVA）结合Bonferroni校正，确定各组织间表达差异的显著性，使用星号标注显著水平，其中****表示P≤0.000 1。

全基因组关联分析所用的303份降香黄檀植株群体栽植于海南大学儋州校区，连续测定了12个月份的株高数据。基于全基因组重测序数据，使用GEMMA软件结合亲缘关系与群体结构，采用混合线性模型（Mixed Linear Model，MLM）对303份样本的株高表型数据进行关联分析，显著性阈值采用Bonferroni校正法确定，设置P value<1×10⁻⁵的SNP（Single Nucleotide Polymorphism，SNP）为显著的SNP。

选取降香黄檀叶片cDNA为模板，利用同源重组方法将扩增片段构建至pCAMBIA1300载体中。将构建好的pCAMBIA1300-DoMYB48质粒转入GV3101菌株中，在LB液体培养基中进行振荡培养，直至OD₆₀₀达到0.5～0.6，离心收集菌体，使用由10 mmol·L⁻¹ MgCl₂、10 mmol·L⁻¹ MES、100 μmol·L⁻¹ 乙酰丁香酮（Acetosyringone，AS）组成的侵染液重悬菌体，室温暗培养2～3 h。将重悬液注射到烟草叶片中，黑暗培养12 h后光照培养24 h，在激光共聚焦显微镜（LSM900）下对荧光信号进行亚细胞定位情况判断。

本研究使用Pfam数据库中MYB基因的HMM模型在降香黄檀基因组中进行筛选，经Pfam及NCBI结构域验证并去除冗余后，共筛选鉴定到180个MYB家族成员（DoMYB1～DoMYB180）。根据R结构域数量分为1R-MYB（38个）、2R-MYB（137个）、3R-MYB（5个）3种类型，未检测到4R-MYB类型成员。染色体定位显示180个DoMYB基因不均匀地分布于全部10条染色体上，不同染色体上的基因数量介于14至29之间，其中2号、3号、8号及9号染色体上存在明显的基因簇（图1）。

Figure 1.  Gene distribution of the DoMYB gene family on D. odorifera chromosomes

理化性质分析表明，DoMYB蛋白相对分子质量范围为8.17～183.08 kDa；氨基酸数目差异显著，最少的有72个氨基酸，最多的有1 683个；等电点范围为4.27～9.92；不稳定系数介于30.86～88.91之间，仅有6个蛋白的不稳定系数小于40，多数DoMYB蛋白为不稳定蛋白；脂肪系数范围为50.00～91.25；除DoMYB130外，其他蛋白的平均亲水性小于0，表明该家族成员多为不稳定亲水性蛋白。3R-MYB亚类的平均氨基酸数目和分子质量均高于1R-MYB和2R-MYB，其余理化性质在各亚类间无显著差异。亚细胞定位预测显示，全部成员表现出细胞核特异性定位的特点。

以180条DoMYB蛋白序列为基础构建系统发育树，根据系统发育关系将DoMYB基因划分为23个亚组（A1～A23）。其中，1R-MYB的38个成员集中分布于A10、A11、A13、A14和A21亚组中，其余分散在A1、A7、A9、A18、A22和A23亚组中；3R-MYB主要分布在A22亚组中，2R-MYB类型的137个成员分布最为广泛，遍布多个亚组。

保守结构域分析显示，所有的DoMYB基因都包含典型的MYB结构域（PLN03091、PLN03212、SANT），其中PLN03091占比最高（97个，53.9%），其次为PLN03212结构域（39个，21.7%）与SANT结构域（18个，10.0%）。其中DoMYB9、DoMYB105及DoMYB122不仅含有典型的MYB结构域，还包含COG5114、RSC8等其他类型的结构域（图2-A）。MEME分析鉴定出20个不同的保守基序（Motif1～Motif20），每个蛋白包含1到10个不等的基序。Motif3（91.7%）、Motif1（85.0%）和Motif5（77.8%）出现频率最高，其余的基序出现频率为3.9%～26.1%。同一亚组内的基因通常具有相似的基序组成，如A5亚组均含Motif1、Motif2、Motif3、Motif5和Motif12，A11亚组中均含Motif8、Motif10、Motif18。此外，亚组中存在不同分支具有不同Motif数目的情况，A22亚组10个基因均含Motif1、Motif2、Motif3、Motif5、Motif12，其中分支1的5个基因（DoMYB8、DoMYB29、DoMYB121、DoMYB141、DoMYB158）相较于同亚组分支2（DoMYB9、DoMYB42、DoMYB48、DoMYB122、DoMYB151）多出Motif8（图2-B）。基因结构解析结果表明，外显子数量范围为1～12个，内含子数量范围为0～11个，其中120个基因含3个外显子，29个基因含有2个外显子，仅有9个基因含有1个外显子，其余的22个基因的外显子数量范围为4～12个。同亚组中亲缘关系紧密的基因，其基因结构模式较为一致，而不同亚组间UTR、CDS和内含子在数量和位置上存在差异（图2-C）。

Figure 2.  Analysis of conserved domains, conserved motifs and gene structure of the DoMYB gene family

基于180个降香黄檀MYB蛋白与168个拟南芥的MYB蛋白序列共同构建系统发育树，被划分为25个亚组（C1～C25）（图3）。1R-MYB成员主要分布在C1～C5亚组，3R-MYB蛋白主要集中于C23亚组，2R-MYB蛋白分布于整个进化树范围内。各亚组成员数量差异显著，C5亚组最多（60个），C13亚组（37个）和C12亚组（35个）次之。根据拟南芥中MYB转录因子的已知功能，可推测降香黄檀同源基因功能。C9亚组中的DoMYB112、DoMYB113、DoMYB148、DoMYB97、DoMYB90、DoMYB27等与拟南芥基因AtMYB123聚类于同一分支中，这些基因可能在根系发育与形态建成发挥作用^[24]；C20亚组中的DoMYB48与AtMYB104聚类于同一分支，属于GAMYB亚类成员，可能在植物的种子萌发、茎的伸长及花药分化过程中发挥重要的作用^[25]。

Figure 3.  Phylogenetic analysis of the MYB gene families in D. odorifera and A. thaliana

对DoMYB基因进行种内共线性分析，共鉴定出56对共线性关系的基因对。主要分布在2号、3号及5号染色体中，其中2R/2R类型基因对46对，1R/1R类型基因对4对，1R/2R类型基因对4对，3R/3R类型基因对2对。计算了56对重复基因对的Ka、Ks以及Ka/Ks比值以评估选择压力，结果显示所有可分析基因对的Ka/Ks比值均未超过1（0.106～0.465），表明这些基因对受纯化选择作用（图4、表2）。

Figure 4.  Intra-species collinearity analysis of DoMYB genes

进一步将降香黄檀与拟南芥、大豆、水稻和毛果杨4个物种进行种间共线性分析，降香黄檀与4个物种均存在共线性基因对，与拟南芥之间鉴定出201对同源基因对，与大豆之间鉴定出699对同源基因对，与水稻之间鉴定出92对同源基因对，与毛果杨之间鉴定到386对同源基因对。进一步分析表明，DoMYB133与其他物种的20个基因存在同源共线性，该基因与拟南芥中5个基因同源，与大豆中8个基因同源，与毛果杨中6个基因同源，与水稻中1个基因同源，是具有最多共线性对的基因。此外，DoMYB19、DoMYB53、DoMYB98等47个基因在拟南芥、大豆、水稻和毛果杨之间都存在共线性关系，表明这些基因在MYB基因家族的长期进化过程中发挥了关键的调控作用（图5）。

Figure 5.  Interspecific syntenic relationships of MYB family genes among D. odorifera, A. thaliana, soybean, rice, and P. trichocarpa

为探究降香黄檀DoMYB基因转录调控潜力，对其上游启动子区域内的顺式作用元件进行预测，共筛选到41种类型各异的顺式调控功能元件，对应的元件总数为4 541个（图6）。将这些功能序列根据其生物学功能归属，划分为4大功能类别：光响应元件（13种）、生长发育相关元件（12种）、植物激素响应元件（10种）以及生物及非生物胁迫响应元件（6种）。光响应元件的占比与数量高于其他类别，179个基因含有光响应元件，其中Box 4和G-box出现频率较高，表明DoMYB基因家族普遍受到光照信号的调控。胁迫响应元件以厌氧诱导元件（ARE）分布最广（136个基因），耐旱元件（MBS）和低温响应元件（LTR）分别分布于约1/3和1/4的成员中。生长发育相关元件如分生组织表达元件（CAT-box）分布于63个基因中，昼夜节律控制元件（Circadian）分布于42个基因中。此外，植物激素响应元件中，脱落酸响应元件（ABRE）分布于136个基因中，茉莉酸甲酯响应元件（TGACG-motif、CGTCA-motif）分布于125个基因中，响应赤霉素元件（TATC-box、GARE-motif、P-box）分布于90个基因中。上述结果表明，DoMYB基因家族的启动子区域富含各类顺式作用元件，表明其成员在多个生物学过程中可能发挥重要的转录调控功能。

Figure 6.  Analysis of the cis-acting elements of the promoter regions of MYB gene family members in D. odorifera

利用降香黄檀3年生植株的韧皮部（B01～B05）、边材（S01～S05）、过渡区（T01～T05）及叶片（L01～L05）组织的转录组测序数据，对DoMYB在各组织样本中的表达模式进行统计解析，综合其相对组织特异性及表达量数据，不同成员呈现显著的组织差异表达（图7）。在叶片中的高表达基因数量最多，占总基因数的52.7%，其中DoMYB56和DoMYB155在叶片中的表达水平显著高于其他组织；韧皮部高表达基因占总基因数的15.6%，包括DoMYB83、DoMYB140、DoMYB159等，其中DoMYB159在韧皮部中表现出明显的组织特异性；边材高表达基因占总基因数的15.6%，包括DoMYB30、DoMYB65、DoMYB79、DoMYB169等；过渡区高表达基因占总基因数的7.8%，包括DoMYB6、DoMYB84、DoMYB89、DoMYB138、DoMYB179等。DoMYB基因在不同组织中发挥重要且独特的作用，可能参与组织特异性的基因调控网络。降香黄檀种内发育树中的A22亚组中包含2R-MYB成员和3R-MYB成员，其两个分支在保守基序和表达模式上存在显著差异。分支1中DoMYB8、DoMYB29、DoMYB121在边材和过渡区中高表达，DoMYB141和fgvbDoMYB158在叶片中高表达；而分支2基因整体主要在叶片中高表达，在其他组织中表达受到抑制，其中DoMYB48在L01与L02样本中高表达。

Figure 7.  Analysis of the expression pattern of the DoMYB gene family

前期课题组通过全基因组关联分析，筛选到与株高显著关联的DoMYB48基因（P=4.151×10⁻⁶）。为进一步研究DoMYB48在不同组织样本中的表达模式，采用RT-qPCR检测该基因在7年生降香黄檀韧皮部、边材、过渡区以及叶片中的表达量，结果显示，DoMYB48在7年生降香黄檀过渡区组织中表达量最高，与叶片组织中的表达量相比呈现出583倍的显著差异（图8）。

Figure 8.  Tissue-specific expression of the DoMYB48 gene in 7-year-old D. odorifera

为了进一步探究DoMYB48基因功能，本研究通过构建重组pCAMBIA1300-eGFP-DoMYB48载体，转化农杆菌后瞬时转化至烟草叶片中，在激光共聚焦显微镜下观察DoMYB48融合蛋白的荧光信号，并以pCAMBIA1300-eGFP空载作为对照。注射空载菌液的烟草叶片中，绿色荧光在细胞核和细胞膜上分布；注射pCAMBIA1300-eGFP-DoMYB48载体菌液的烟草叶片中，绿色荧光分布在细胞核中，与细胞核标记蛋白的红色荧光重叠，表明DoMYB48蛋白定位于细胞核，符合转录因子的调控功能特征（图9）。

Figure 9.  Subcellular localization analysis of DoMYB48

MYB转录因子广泛参与各类生长发育进程、次生代谢、逆境响应等关键生物学过程^[26]。本研究在降香黄檀中共鉴定到180个DoMYB成员，分为1R-MYB（38个）、2R-MYB（137个）和3R-MYB（5个）三种类型。与拟南芥（198个）^[27]、大豆（244个）^[28]、毛果杨（192个）^[29]等物种数量相当，反映了MYB家族在植物进化中的保守性；与小麦（Triticum aestivum）（719个）^[30]、棉花（Gossypium hirsutum）（812个）^[31]等相比，MYB基因家族的数量差异反映在进化过程中全基因组多倍体化事件驱动的MYB家族扩张变化规律。DoMYB基因在10条染色体上呈不均匀分布，部分基因呈簇状分布，这种分布与串联重复事件相关，是基因家族扩张的重要机制^[32]。系统发育分析将180个DoMYB成员划分为23个亚组，与拟南芥MYB家族的分类系统具有良好的对应关系。基因组内共鉴定出56对共线性基因对，Ka/Ks比值均小于1，重复基因经历了强烈的纯化选择后得以保留^[33]。种间共线性分析显示，降香黄檀与大豆的共线性对远多于与水稻的共线性对，这与豆科内部的近缘关系一致。47个DoMYB基因在拟南芥、大豆、水稻和毛果杨中均存在同源基因，可能具有潜在的跨物种功能保守性，在MYB家族进化中具有重要作用。

理化性质分析表明，家族成员多为不稳定亲水性蛋白，其氨基酸长度差异显著，该家族可能参与从基础代谢到发育等广泛与复杂的调控^[34]。所有成员均预测定位于细胞核，这与MYB转录因子在核内调控基因表达的功能一致。对保守基序进行分析显示，20个保守基序在不同基因间呈现不均匀分布，具有不同保守基序的基因可能参与特定的生理途径或胁迫响应。A22亚组内出现基序分化：分支1含有Motif8，而分支2缺失该基序，基序的丢失可能导致DNA结合特异性或蛋白互作能力的改变，进而驱动功能分化^[35]。基因结构分析显示，多数基因的结构为3个外显子和2个内含子的模式，这与油菜（Brassica napus）中的研究一致^[36]。启动子区域顺式作用调控元件分析显示，光响应元件最为丰富，表明MYB基因在光信号响应及光形态建成中扮演重要角色^[37]。

组织表达分析显示，大多数DoMYB基因在叶片中的表达量较高，这与启动子区域富含光响应元件结果相符，也与叶片作为光合作用主要场所的生理功能相一致。在叶片中高表达的DoMYB155位于C9分支，与AtMYB123聚类，可能在花青素合成与种子发育中发挥关键作用^[24]。在韧皮部中高表达的DoMYB159位于C13分支，与HAG亚家族的AtMYB28和AtMYB29聚类，其在拟南芥中参与硫代葡萄糖苷的生物合成调控，表明DoMYB159可能在韧皮部防御反应或次生代谢物合成中发挥重要作用^[38]。在边材中高表达的DoMYB30、DoMYB65、DoMYB79和DoMYB169位于C15分支，与AtMYB26、AtMYB46、AtMYB83、AtMYB103等次生细胞壁及木质素生物合成重要基因聚类，可能参与边材的次生细胞壁形成和木质化过程^[39-40]。在过渡区中高表达的DoMYB179与AtMYB46、AtMYB83、AtMYB103等次生细胞壁相关基因聚类于C15分支，可能参与边材向心材过渡区中细胞壁组分的调控与代谢物积累过程^[41]。不同DoMYB基因分别调控相应组织中的特异基因表达网络，为理解MYB转录因子在植物生长发育中的复杂调控机制提供了重要线索。

基于课题组前期研究，筛选到与株高显著相关的2R-MYB成员DoMYB48，该基因与拟南芥C20亚组中AtMYB104聚类，属于与赤霉素调控相关的GAMYB分支，参与花药发育及GA信号传导，调节茎伸长及植株开花^[42-44]。组织表达分析显示，DoMYB48在3年生降香黄檀植株中叶片表达量最高，可能主要参与营养生长阶段的叶片发育及茎伸长过程；在7年生植株中过渡区表达量最高，表明其功能可能随发育阶段转换，从早期调控叶片发育转向参与心材形成的次生代谢调控，这种转变可能与降香黄檀心材形成的树龄相关，7年生植株进入心材发育阶段，过渡区成为代谢活跃区^[45]。在其他木本植物中，MYB转录因子随着发育阶段转换的现象已有报道；在杉木中，基因在幼龄茎部表达较高，但随着茎的成熟表达量逐渐下降^[46]。DoMYB48基因的表达模式随着植株的发育过程产生了变化。

综上所述，在降香黄檀基因组中共鉴定到180个DoMYB基因家族成员，在保守基序、基因结构及理化性质方面表现出明显的特异性与多样性。启动子区域内光响应元件最为丰富。组织表达分析显示DoMYB基因呈现组织特异表达模式，其中DoMYB48在降香黄檀的不同发育阶段具有不同的差异表达模式。本研究为解析MYB转录因子在降香黄檀生长发育中调控机制及后续功能验证提供了理论基础。