基于全长转录组测序的掌叶大黄bZIP基因家族鉴定分析

唐 璟，杜 桥，胡晓晨，李依民，高 静，彭 亮，张 岗

（1.陕西中医药大学药学院 / 陕西省秦岭中草药应用开发工程技术研究中心, 陕西 西安 712046；2.陕西中医药大学省部共建特色秦药资源研究开发国家重点实验室(培育), 陕西 咸阳 712083；3.陕西中医药大学陕西省中医药管理局“秦药”研发重点实验室, 陕西 西安 712046）

bZIP (basic leucine zipper)转录因子是真核生物中数量最大、最保守的基因家族之一，该家族成员均含有由一个碱性区域 (basic region) 和一个亮氨酸拉链结构域 (leucine zipper) 组成的特征保守结构域[1]。碱性区域含有一个核定位信号和一个保守的N-X7-R/K 结构，可特异性结合DNA；而亮氨酸拉链区域为L-X6-L-X6-L 结构，具有二聚化能力，该区域形成二聚体后依赖于碱性区域结合在双链DNA 分子上发挥作用[2]。bZIP 转录因子可特异性结合DNA启动子区的ACGT 核心顺式作用元件，如ABRE(ABA-responsive element)、A-box (TACGTA)、C-box(GACGTC)、G-box (CACGTG)，从而调节下游基因的表达[3-4]。

许多bZIP 转录因子在植物体内的生物功能已被阐明，广泛参与调控植物体生长发育、次生代谢物合成及逆境胁迫等生物学过程。在拟南芥(Arabidopsis thaliana) 中，bZIP17与bZIP28一起介导细胞生长的多个基因的非诱导表达，从而调控根的伸 长[5]。ABA 反 应 元 件 结 合 因 子 (ABA-responsive element binding factors)ABF2、ABF3、ABF4通过转录激活叶绿素分解代谢基因和衰老相关基因，促进脱落酸 (abscisic acid, ABA) 介导的叶绿素降解和叶片衰老[6]。TabZIP74正调控小麦 (Triticum aestivum) 条锈病抗性[7]。此外，bZIP 转录因子参与激素信号转导，激活某些次生代谢物的生物合成通路基因表达，影响次生代谢物的合成。葡萄 (Vitis vinifera)VvibZIPC22可激活类黄酮合成途径关键基因的启动子进而指导黄酮醇生物合成[8]。DkbZIP5通过结合DkMYB4启动子区域的ABRE 元件，并以依赖ABA 的方式，调控参与原花青素生物合成的基因DkMYB4的 表 达[9]。 TGACG 结 合 因 子 (TGACG motif-binding factors)TGA6通 过 响 应 水 杨 酸 信 号 调控黄花蒿 (Artemisia annua) 中青蒿素的生物合成[10]。AabZIP1在ABA 诱导下直接结合AaMYC2的启动子而上调其表达，从而正向调控青蒿素的生物合成[11]。同时，药用植物bZIP基因通过介导MeJA 信号转导通路参与次生代谢物合成。雷公藤 (Tripterygium wilfordii)TwTGA1受 茉 莉 酸 甲 酯 (methyl jasmonate，MeJA) 诱导，能调节长春花生物碱合成[12]。在MeJA 诱导下，bZIP17和bZIP60参与调控紫花苜蓿 (Medicago sativa) 三萜皂苷生物合成[13]。

蓼科大黄属多年生高大草本植物掌叶大黄(Rheum palmatum) 是我国传统大宗药材大黄的基源之一，因其喜凉爽湿润气候，故多分布于中国甘肃、四川、青海、陕西等高海拔区域，今多为栽培品[14]。大黄以干燥根及根茎入药，性味苦、寒，具有泻下攻积、清热泻火、凉血解毒、逐瘀通经、利湿退黄等功效[15]。大黄含有蒽醌类、黄酮类等化学成分，其中蒽醌类和蒽酮类为大黄的特征性成分，也是大黄发挥药效的主要活性成分[16]。次生代谢产物生物合成途径复杂，通过转录因子调控其中的一系列关键酶基因的表达，进而影响药材品质形成[17-18]。鉴于bZIP 转录因子对植物次生代谢产物合成的调控作用，掌叶大黄bZIP 转录因子的鉴定和功能研究尚未报道。为此，本研究基于三代全长转录组测序挖掘掌叶大黄bZIP 转录因子家族，对该家族成员进行特征分析，同时结合RNA-seq 数据分析RpbZIPs基因在不同器官及MeJA 处理下的表达模式，为进一步研究大黄bZIP 转录因子的生物学功能奠定基础。

1 材料与方法

1.1 主要仪器和试剂

MeJA、无 水 乙 醇；RN38-EASYspin Plus 植 物RNA 提 取 试 剂 盒 (Aidlab 公 司)； PrimeScript™RT Master Mix 反 转 录 试 剂 盒、TB Green®Premix ExTaq™ Ⅱ (TliRNaseH Plus) (TaKaRa 公司)。K5800 自动检测超微量分光光度计 (凯奥公司)、StepOnePlus™Real-Time PCR (qPCR) 仪 (美 国Applied Biosystems公司)。

1.2 材料

掌叶大黄一年生植株和成熟种子于2021 年8 月采集于甘肃中医药大学和政药用植物园，经陕西中医药大学胡本祥教授鉴定。

取一年生植株3 株，将叶、根、根茎3 个部位各等量混合后，由广州基迪奥生物科技有限公司测序分析。在每个装有200 g 泥炭土 (klasmann) 的塑料花盆中 (黑色，直径9 cm、高12.5 cm)点播4 颗大小均匀、饱满的掌叶大黄种子，培养条件光周期为光照16 h、黑暗8 h，光照强度9 000 lx，温度(23 ± 2) ℃。首次浇水200 mL，之后每3 d 浇水50 mL。选取完整均匀的一月龄幼苗，激素处理组对整株喷施200 μmol·L-1MeJA，对照组喷施溶剂，所有样品进行生物学重复3 次，分别于处理后12 h 取样，将根和叶等量混合于液氮中速冻，进行RNA-Seq 分析。

选取长势均匀、完整的一月龄幼苗，处理组喷施200 μmol·L-1MeJA，对照组喷施溶剂(将助溶剂95%乙醇与试验组相同倍数进行稀释)；两组均设置0 h 为空白对照，Mock 为溶剂对照，所有样品重复3 次，分别于处理后3、6、12 和24 h 进行取样，于液氮速冻后置-80 ℃冰箱保存，供后续qRT-PCR 验证分析。

1.3 全长转录组测序及RpbZIP 鉴定

掌叶大黄全长转录组测序数据为课题组前期所构建[19]。使用BlastX (http://www.ncbi.nlm.nih.gov/BLAST/) 将全长转录组测序所得完整的isoform 在Nr (http://www.ncbi.nlm.nih.gov)、SwissProt (http://www.expasy.ch/sprot) 和KEGG (http://www.genome.jp/kegg)数据库中进行功能注释分析。

将预测的蛋白序列与Plant TFDB (http://planttfdb.gao-lab.org/) 进行 hmmscan 比对，将预测到的TF 进行归类后，从中筛选bZIP 转录因子家族基因。通过BlastX 比对、ORF Finder 分析筛选出具有完整ORF的bZIP全长基因。利用NCBI-CD search 和ExPASy(https://prosite.expasy.org/prosite.html) 验 证 编 码 蛋 白的保守结构域。

1.4 蛋白序列特征分析

利用ProtParam (https://web.expasy.org/protparam/)和SOPMA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_autom at.pl?page=npsa_sopma.html) 对 掌 叶 大 黄bZIP 转 录因子蛋白质理化性质进行分析；用Cell-PLoc 2.0(http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) 进行亚细胞定位；利用软件Jalview (https://www.jalview.org/download/)的MAFFT 功能进行多序列比对；MEME(https://meme-suite.org/meme/tools/streme) 进行Motif分析，并用TBtools 进行可视化。

1.5 系统进化分析

以已鉴定完成的拟南芥bZIP 氨基酸序列为参照，使用MEGA 7 将掌叶大黄与拟南芥的bZIP 蛋白序列进行系统进化分析 (采用邻接法，Bootstrap 为1 000 次)。

1.6 基因表达模式分析

基于掌叶大黄不同部位的RNA-Seq 转录组测序数据，分析bZIP 转录因子基因表达差异 (以FPKM 值计算)；另从MeJA 处理的掌叶大黄转录组数据中筛选显著差异表达的RpbZIP基因 (FDR ＜0.05 且|log2FC| ＞ 1)。二者均用TBtools 进行均一化和可视化绘图。

1.7 qPCR 验证

基于RNA-Seq 转录组测序数据选择MeJA 处理后表达显著上调或下调的RpbZIP基因进行qPCR验证，以β-actin[20]为内参基因。利用NCBI-Primer(https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi )设计引物 (表1)，由北京奥科合成引物。用RN38-EASYspin Plus 植物RNA 快速提取试剂盒提取植物总RNA，1.0%琼脂糖凝胶电泳和K5800 超微量分光光度计用于检测RNA 的质量。用反转录试剂盒(PrimeScript™ RT Master Mix) 合 成cDNA 第 一 链。qPCR 反应按照TB Green®Premix ExTaq™ (TaKaRa,RR420A, 中国) 试剂盒说明书配置10 μL 体系。反应程序：预变性95 ℃ 30 s，变性 95 ℃ 5 s，退火 60 ℃30 s，延伸 60 ℃ 34 s, 40 个循环。技术重复和生物重复各3 次，应用2-ΔΔCt法计算基因的相对表达量。

表1 掌叶大黄bZIP 基因qPCR 引物Table 1 The qPCR primers for bZIPs genes in Rheum palmatum

1.8 数据处理

利用SPSS 26 软件对掌叶大黄bZIP基因响应外源激素的表达量进行方差分析和显著性检验，并使用Graphpad 8.0 绘图。P＜ 0.05 代表有显著性差异。

2 结果与分析

2.1 RpbZIPs 基因鉴定及蛋白理化性质分析

从掌叶大黄全长转录组中鉴定到63 个具有完整结构域的全长RpbZIPs，按照RNA-Seq 中的isoform顺序将其重命名为RpbZIP1～RpbZIP63。RpbZIPs蛋白理化性质分析发现，RpbZIPs基因编码蛋白的氨基酸残基数目为144 (RpbZIP54)～787 (RpbZIP1)aa (表2)，分子量16 457.34～85 075.16 Da，等电点为5.04 (RpbZIP55)～10.10 (RpbZIP59)，说明主 要为中性和碱性蛋白；RpbZIP10和RpbZIP34所编码的蛋白不稳定系数分别为27.80 和30.23，为稳定蛋白，其余蛋白不稳定系数均大于40，为不稳定蛋白；脂肪系数为50.84 (RpbZIP7)～87.62 (RpbZIP19)；总平均亲水性均为负值，均为亲水性蛋白。除了RpbZIP60和RpbZIP63没有β 转角，其余61 个RpbZIPs均由α 螺旋、延伸链、β 转角和无规卷曲结构组成，且以α 螺旋和无规则卷曲为主。

2.2 RpbZIPs 转录因子系统发育分析

利用MEGA 7.0 构建RpbZIPs 与拟南芥AtbZIPs转录因子的进化关系 (图1) 。根据拟南芥的分组规则将所有bZIPs 分为13 个亚族 (A～K 和M、S)，除了J 和M 亚族无掌叶大黄bZIP 成员，其他11 个亚族都同时含有拟南芥和掌叶大黄bZIP 成员。其中，A 组成员最多，包含RpbZIP 成员12 个，AtbZIP 成员13个；I 组所含RpbZIP 成员最多，为14 个，而拟南芥有10 个成员，掌叶大黄的I 家族成员可能发生复制事件；M 和J 组均只包含1 个AtbZIP 成员。

图1 掌叶大黄与拟南芥bZIP 成员系统发育关系Figure 1 Phylogenetic relationships of Rheum palmatum and Arabidopsis bZIP members

2.3 RpbZIPs 蛋白结构分析

RpbZIPs 多序列比对结果表明，63 个RpbZIPs蛋白都包含典型的bZIP 家族结构域 (图2) ，即碱性区域和亮氨酸拉链结构域。该碱性结构域包含一个保守的N-X7-R/K 结构，其次是两个保守的碱性氨基酸精氨酸 (R) 和赖氨酸 (K)；亮氨酸拉链结构域的特征为L-X6-L-X6-L，碱性区域和第一个亮氨酸相对保守。此外，发现保守区域部分氨基酸被疏水残基取代，其中RpbZIP8 碱性区域的R / K 被一个异亮氨酸 (I) 取代；亮氨酸拉链的部分亮氨酸被异亮氨酸、缬氨酸 (M)、甲硫氨酸 (F) 等疏水氨基酸取代。

图2 RpbZIPs 多序列比对Figure 2 Sequence alignment of RpbZIPs

使用MEME 对63 个RpbZIPs 蛋白进行保守基序分析 (图3)，共鉴定了12 个Motif，所有RpbZIP 蛋白均含有典型碱性亮氨酸拉链基序Motif1，Motif2只存在于I 和E 亚族，Motif3 只存在于S 和C 亚族，Motif4 为I 亚族特有基序，Motif6、Motif8、Motif12仅于A 亚族部分成员中发现。同一亚族内RpbZIPs所含Motif 的种类和总体分布相似，同一亚族内的RpbZIPs 可能有相似的功能。

2.4 RpbZIPs 基因表达模式分析

基于转录组数据，对63 个RpbZIPs在掌叶大黄叶、根及根茎3 个部位的基因表达进行热图聚类分析。RpbZIPs基因表达具有组织特异性，23 个RpbZIP基因主要在叶中表达 (图4A)；39 个RpbZIP基因主要在根或根茎中表达；RpbZIP49在3 个部位中均无表达。其中有部分处于同一亚族的RpbZIP基因的表达模式相似，可能在相同的部位共同发挥功能，如E 亚族的RpbZIP29、RpbZIP36、RpbZIP40、RpbZIP42、RpbZIP48和H 亚族的RpbZIP60、RpbZIP63等。

图4 掌叶大黄不同组织及MeJA 处理12 h 后RpbZIP 基因的表达热图Figure 4 Expression heat map of RpbZIP genes in different tissues and RpbZIP genes at 12 h after MeJA treatment

基于MeJA 处理12 h 的大黄幼苗的转录组数据，筛选出13 个显著差异表达的RpbZIPs基因。热图聚类分析结果表明，7 个RpbZIP基因 (RpbZIP14、RpbZIP16、RpbZIP17、RpbZIP30、RpbZIP41、RpbZIP43、RpbZIP61) 在MeJA 处理12 h 后基因表达上调 (图4B)，6 个RpbZIP基 因 (RpbZIP19、RpbZIP23、RpbZIP26、RpbZIP36、RpbZIP42、RpbZIP59) 表达下调。同时发现，RpbZIP16、RpbZIP17、RpbZIP61等基因在大黄药用部位根和根茎中的表达量较高且响应MeJA。

2.5 RpbZIPs 候选基因的表达验证

以MeJA 处理0 、3 、6、12、24 h 的掌叶大黄幼苗为材料，以MeJA 差异表达候选的6 个基因为候选基因，进行RpbZIP基因表达qPCR 验证 (图5)。以0 h (CK) 作 为 空 白 对 照，Mock 作 为 溶 剂 对 照，RpbZIP14、RpbZIP16、RpbZIP17、RpbZIP61在 处 理12 h 的表达量高于溶剂对照组，RpbZIP26、RpbZIP42在MeJA 处理12 h 的表达量低于溶剂对照组，与转录组测序基本一致。

图5 MeJA 处理下RpbZIP 基因的相对表达量Figure 5 Relative expression levels of RpbZIP genes under MeJA treatment

RpbZIP14、RpbZIP16、RpbZIP61在MeJA 处理3 h的表达量均达到峰值，RpbZIP14在6 h 开始表达量呈下降再上升趋势；RpbZIP16在3 h 表达量达到峰值为CK 的71.23 倍，之后呈下降趋势；RpbZIP17的表达量在12 h 内持续上调至峰值为CK 的68.75倍；RpbZIP26在MeJA 处理3～24 h 的表达量均低于CK 的0.5 倍；RpbZIP42在24 h 内(除6 h 外)表达量持续下调。

3 讨论

bZIP 转录因子通常以二聚化形式依赖碱性区域特异性结合DNA 启动子区，从而参与基因的表达调控，控制植物体内多种生命过程[21-22]。随着转录组与基因组测序技术的发展，bZIP 在拟南芥[23]、水稻 (Oryza sativa)[24]、玉米 (Zea mays)[25]等植物中已被广泛研究。本研究基于掌叶大黄全长转录组数据，鉴定出63 个RpbZIPs 转录因子基因。其数目少于药用植物丹参 (Salvia miltiorrhiza) (70)[18]、紫花苜蓿 (138)[26]、薏 苡 (Coix lacryma-jobi) (81)[27]、野 菊(Chrysanthemum indicum) (75)[28]等，多于大麻 (Cannabis sativa) (55)bZIP数目[29]，说明bZIP基因家族数目在不同物种中变化较大。进一步结合掌叶大黄根、根茎及叶片3 个部位及MeJA 处理的比较转录组数据分析，富集到药用部位中大量表达且同时受MeJA诱导的候选基因，并验证了6 个基因与转录组测序结果基本一致。本研究为后续RpbZIP 功能研究及其调控机制奠定了基础。

RpbZIP 蛋白主要为不稳定亲水蛋白，等电点5.04～10.10，说明RpbZIP 转录因子在多种微环境中发挥作用，二级结构以α-螺旋和无规则卷曲为主。系统进化表明，除J 和M 亚族无RpbZIP 成员外，RpbZIP 家族成员分布在11 个亚族内。多序列比对发现RpbZIP 的碱性结构区域最保守，而亮氨酸区域含有大量的疏水残基，符合bZIP 的结构特征[22]。除了共有保守结构域，RpbZIP 不同亚族之间还存在其他保守位点，其中RpbZIP8 碱性区域的R / K 被一个异亮氨酸 (I) 取代，该特征与黄花蒿中AabZIP78的结构变化相一致[30]；E 亚族RpbZIPs 在亮氨酸拉链区域有一个保守的脯氨酸残基，可干扰同源二聚体的形成[31]；I 亚族成员在R / K 保守位置均为K，该特征可能与维管束发育功能相关[22]。MEME 分析显示，同亚族RpbZIP 成员拥有相似Motif，结合多序列比对，存在于所有RpbZIP 中的Motif 1 很可能代表bZIP 保守结构域，而其他Motif 可能影响基因家族的功能多样性[32]。

bZIP 各个亚族广泛参与植物生长发育过程，A亚族主要参与ABA 响应和逆境胁迫的转录调控[33]；D 亚族 (TGA 因子) 主要参与防御反应及根的生长发育[34-35]；F 亚族参与植物锌防御反应和盐胁迫响应[36]；H 亚族主要由HYH 和HY5 转录因子组成，参与次生代谢产物生物合成与光信号的转导等[37]；而S 亚族的生物功能较为广泛，参与植物生殖发育，还可响应机械刺激、干旱和低温等[22]。H 亚族RpbZIP成员均在叶中表达，根及根茎中几乎无表达，推测这些基因主要参与叶中次生代谢产物生物合成与光信号的转导。S1 组AtbZIP11控制生长素依赖的主根生长[38]；AtbZIP1及AtbZIP53在根中控制盐胁迫下的代谢重编程[39]，而RpbZIP61在掌叶大黄根及根茎中高表达，且受MeJA 诱导，可能参与根及根茎的发育和代谢。

大量研究表明，bZIP 转录因子可响应各种激素，参与植物生物合成和抗逆相关基因的调控。茉莉酸甲酯作为一种重要的植物激素，在植物防御反应及次生代谢产物生物合成中起着重要的作用[40]，常被用于药用植物次生代谢物合成关键转录调控基因挖掘研究。如穿心莲 (Andrographis paniculata)中鉴定发现了7 个ApbZIPs在叶片中表达最高，且MeJA 处理后显著上调，可能参与调控穿心莲内酯的生物合成[41]。本研究利用MeJA 处理12 h 掌叶大黄幼苗转录组数据筛选到13 个显著差异表达的RpbZIPs，结合不同组织部位表达分析，发现RpbZIP16、RpbZIP17、RpbZIP26、RpbZIP61等在掌叶大黄根及根茎中大量表达且响应MeJA 处理。RpbZIP16、RpbZIP17属 于 CPRF1 (Common Plant Regulatory Factor 1) 因 子，受MeJA 诱 导。黄 花 蒿CPRF1 因子AabZIP33和AabZIP38调控蓝光介导青蒿素的生物合成[29]；欧芹 (Petroselinum crispum) CPRF1抑制查尔酮合酶 (CHS) 基因的启动子活性，从而下调黄酮类化合物的积累[42]。RpbZIP26属于TGA 因子，MeJA 处理抑制其表达。大多数TGA 因子含有MeJA 反应元件TGACG-基序，在MeJA 处理下主要参与次生代谢产物生物合成及防御相关基因的表达调控[43-44]。因此，响应MeJA 的RpbZIP16、RpbZIP17、RpbZIP26等可能参与大黄次生代谢产物生物合成的转录调控和抗逆反应，其具体功能还有待于后续深入研究。

4 结论

本研究在掌叶大黄全长转录组数据中鉴定到63 个bZIP 转录因子基因家族成员，并分析其编码蛋白理化性质、多序列比对、系统进化等特征。通过比对掌叶大黄不同组织及MeJA 处理的转录组数据，筛选出掌叶大黄药用部位中高表达且可能响应MeJA 的候选基因，并对其基因表达模式进行qPCR验证，进一步明确RpbZIPs对MeJA 的响应特性，为后续RpbZIP 功能研究及大黄次生代谢转录调控机制研究奠定基础。