干旱胁迫和复水处理后梭梭转录因子的转录组分析

何江峰，王力伟，2，房永雨，王蕴华，王 朝，刘红葵

(1.内蒙古自治区农牧业科学院 生物技术研究中心，内蒙古 呼和浩特 010031；2.内蒙古大学 生命科学学院，内蒙古 呼和浩特 010020；3.内蒙古农业大学 生命科学学院，内蒙古 呼和浩特 010010)

梭梭(Haloxylonammodendron)属黎科(Chenopodiaceae)，为超旱生小乔木，呈高大灌木状，植株高度一般在1～4 m，个别植株可高达8～10 m，亦称盐木、琐琐树、梭梭柴。梭梭是肉苁蓉(CistanchedeserticolaY.C.Ma)的寄主，肉苁蓉是名贵的中药材，在医学上具有抗肿瘤、保肝、补肾、通便、抗老年痴呆、抗辐射等独特功能，其寄生在梭梭的根端[1-2]。梭梭长期生长在荒漠环境中，形成了抗旱、抗风沙、耐高温、耐寒、耐盐碱和耐贫瘠等特性。因此，梭梭是荒漠化地区主要的防风固沙树种，也可作为骆驼、羊适口性饲料，是荒漠草场中主要的木本饲料树种之一[3]。梭梭具有重要的经济价值、生态价值和饲用价值，是我国的二级保护植物，国内外许多学者已从梭梭抗逆的形态学、解剖学、生理学、生态学及分子生物学等方面进行了相关研究。

梭梭的生存环境极其严酷，其主要分布地区常年干旱少雨，年降水量不足50 mm，相对湿度低于40%，蒸发量高于降水量的30～100倍，生存环境的极端温度最高可达40 ℃以上，最低为-40 ℃以下，地表温度可达70 ℃以上[4]。为了适应极端环境，梭梭逐步形成了自身的极端抗旱机制，但梭梭的抗旱机理、抗旱遗传机制和抗旱相关组学的研究相对较少，前人研究发现，部分转录因子及与干旱胁迫相关，例如，AP2-EREBP转录因子(TINY、CBF1、Ptis、AtEBP、DREB1、DREB2)主要调节植物对低温、干旱及高盐[5]等的分子应答反应；bZIP转录因子在拟南芥基因组、豆科植物基因组、水稻基因组中参与对各种逆境胁迫的响应[6-7]；拟南芥中3个NAC基因(ANAC019、ANAC055和ANAC072)的过表达能够显著提高植株对干旱的抗性[8]；许多NAC基因过表达能够显著提高水稻植株的胁迫抗性，却不会减缓植株的正常生长[9-10]；bHLH基因也参与水稻应对干旱胁迫的分子调控过程[11]。在改良植物的抗逆分子育种中，导入或改良转录因子较导入或改良个别功能基因在增强植物抗逆性方面更为有效。因此，本研究主要对干旱胁迫和复水处理下梭梭的转录因子进行分析，挖掘与干旱相关的转录因子，为阐明梭梭抗旱机理奠定理论基础，保护及开发我国抗逆植物基因资源。

1 材料和方法

1.1 材料及处理条件

梭梭种子收集于内蒙古自治区阿拉善左旗。将梭梭种植在塑料花盆中(直径：15 cm，高：18 cm)，盆内土壤均匀，湿度及质量均一致。梭梭幼苗在植物培养室种植，室内温度为24/20 ℃(16 h光照 / 8 h 暗处理)，相对湿度为30%，光子通量密度为350 μmol/(m2·s)。本研究随机单株种植25盆梭梭幼苗，包括1个正常浇水组(对照组，CK)、干旱胁迫组(土壤相对含水量为30%～45%，LMS；土壤相对含水量为20%～25%，MWS；土壤相对含水量为10%～15%，SWS)和旱后复水处理组(土壤相对含水量为在重度水分胁迫组上增加10%～15%，RSWS)。干旱胁迫组、复水处理组和对照组的种植均设置5个生物学重复，每天在固定时间测质量、控水、补水并记录待幼苗生长90 d之后，分别从CK、LWS、MWS、SWS 以及重度胁迫90 d之后复水处理7 d的RSWS的植株上采集等量的根、茎和叶混合，液氮速冻，置于-80 ℃冰箱中储存备用，提取各处理的总RNA。

1.2 总 RNA 提取

材料总RNA的提取采用 RNAprep Pure Plant Plus KIT(DP441，TIANGEN北京生化科技有限公司)试剂盒。用 DNase酶(2212，大连宝生物工程有限公司)对样品中的 gDNA 进行消化处理。RNA完整性和纯度质量经1.2%非变性琼脂糖凝胶电泳和 Agilent 2100 分析系统检测。

1.3 RNA测序和分析方法

总RNA提取按照文献[12]的方法，利用 Solexa/Illumina进行测序及深圳华大基因公司的Illumina系统软件进行数据分析。RPKM值用于基因表达水平的分析，利用Cluster软件[13]，以欧氏距离为距离矩阵计算公式，对基因和试验条件同时进行等级聚类分析，聚类结果用Java Treeview显示，利用WEGO软件对差异基因做GO功能分类统计[14]，pathway显著性富集分析以KEGG pathway为单位，应用超几何检验，找出与整个基因组相比后差异表达基因中显著性富集的pathway基因和基因组的代谢途径[15]。

1.4 de novo 拼接及序列分析

用Trinity[16]做转录组 de novo 组装，并通过 reads overlap 关系得到 Contig。然后，将 clean reads 比对回Contig，通过paired-end reads 能确定来自同一转录本的不同 Contig 以及这些 Contig 之间的距离，Trinity 将 Contig 连在一起，最后得到两端不能再延长的序列，即为 Unigene。组装得到的 Unigene，首先使用 Tgicl 将其去冗余并进一步拼接，然后再对这些序列进行同源转录本聚类，最终得到非冗余 Unigene。用 Tophat v2.0.9(Broad Institute)将拼接得到的非冗余 Unigene 与已公布的同种属基因组信息和转录组数据进行比对、注释，用已获取各测序样本中包含的转录因子编码基因信息。

1.5 差异表达基因的筛选

样本中 Unigene 的表达量以公式化校正的 FPKM(Fragments per kb per million fragments)[17]值衡量。2个样本间差异表达基因(Differential expression genes，DEGs)的筛选及分析使用DEGSeq R package软件，以FDR(False discovery rate)≤0.001，且Log2(fold change，FC)≥1作为筛选阈值进行差异表达基因的筛选。

1.6 差异表达基因 Real-time PCR 验证

从测序结果中随机挑选 9 个转录因子编码基因，用荧光实时定量 PCR 法验证转录组测序数据的有效性和准确性。引物设计用 Primer 5 软件，序列见表1。以18S为内参基因，用2-ΔΔCt计算，并校正目标基因在不同处理下的表达量，每个基因3次生物学重复。10 μL反应体系中包括5 μL SYBR Premix Ex Taq(TaKaRa，TaKaRa Clontech，Japan)，5 μmol/L上游引物，5 μmol/L下游引物和20 ng 模板cDNA。所有反应体系加样都在96孔板进行，并在Lightcycle-480(罗氏)荧光定量PCR仪上完成，反应程序：95 ℃ 5 min；后按95 ℃ 15 s，55 ℃ 2 min，72 ℃ 1 min，共40个循环；72 ℃延伸2 min。18S基因为扩增体系的内参基因，扩增产物的表达水平按 2-ΔΔCt方法计算[18]。每个样品都有3次技术重复，结果都以平均值±标准误(s)(n= 3)计算。

表1 实时定量PCR引物

2 结果与分析

2.1 转录组数据统计分析

以梭梭幼苗为材料，对CK、LWS、MWS、SWS和RSWS处理组进行转录组测序和表达谱分析。对测序数据进行拼接、比对和整理，经de novo 拼接、去冗余、同源转录本聚类后，最终在CK、LWS、MWS、SWS和RSWS 5个组中得到74 641条非冗余Unigene，平均长度780 bp，N50长度1 537 bp。 将CK、LWS、MWS、SWS和RSWS混合样本拼接得到的非冗余Unigene进行比对、注释，共检测到56个转录因子家族，1 307个转录因子编码基因；在CK、LWS、MWS、SWS和RSWS组分别检索到 934，978，994，1 003，991 条转录因子编码基因。LWS、MWS、SWS分别与CK组相比，RSWS与SWS相比，差异表达转录因子均以下调表达方式为主，CK-vs-SWS组中下调转录因子数目最多，达112条，约占该组差异表达转录因子编码基因总数的59.3%(表2)。

表2 差异表达转录因子

2.2 差异表达转录因子家族统计分析

获得4个对比组转录因子基本数量信息后，进一步统计分析转录因子家族及其包含成员的数量关系。LWS、MWS、SWS分别与CK组对比，RSWS与SWS相比，表3分析表明，其中，AP2-EREBP、WRKY、bHLH、MYB、NAC和ABI3VP1 家族富集的DEGs 超过该组差异表达基因总数的50%以上，其中AP2-EREBP、WRKY、MYB和bHLH家族转录因子较多。由图1-A-C结果显示，AP2-EREBP家族包含的下调差异表达基因，分别为6个(CK-vs-LWS)、17个(CK-vs-MWS)和 16个(CK-vs-SWS)DEGs，约占各组DEGs总数的75%，81%，59%；WRKY家族包含的下调差异表达基因，分别为6个(CK-vs-LWS)、12个(CKvs-MWS)和 2个(CK-vs-SWS)DEGs，约占各组DEGs总数的100%，100%，33%；bHLH家族包含的下调差异表达基因，分别为5个(CK-vs-LWS)、4个(CK-vs-MWS)和 16 个(CK-vs-SWS)DEGs，约占各组DEGs总数的100%，57%，84%；MYB家族包含的下调差异表达基因，分别为5个(CK-vs-LWS)、9个(CK-vs-MWS)和 10 个(CK-vs-SWS)，约占各组DEGs总数的100%，64%，45%。图1-D结果显示，AP2-EREBP家族、NAC家族和WRKY家族包含的上调差异表达基因分别占各自DEGs总数的67%，86%，75%。以上表达规律显示，MYB和WRKY家族基因在LWS组中没有上调表达，但在MWS和SWS中其上调表达数逐渐增加。综合说明，RSWS组中偏向通过增加 AP2-EREBP、NAC家族和WRKY家族编码基因的上调表达数来实现梭梭幼苗对水分刺激的应答调节；而干旱组(LWS、 MWS和 SWS)更偏向于下调AP2-EREBP和bHLH 的DEGs数量来实现其对干旱胁迫的响应。

表3 干旱和复水组各转录因子家族中的基因表达数

2.3 差异表达转录因子表达量分析

从4个对比组中挑选了表达量(FC)变化幅度最大的10个上调和下调基因，并对这些基因的表达规律进行统计分析。结果表明：差异表达基因Tify(CL5087)和NAC(Unigene28108)在干旱胁迫组(LWS、MWS和SWS)和复水组(RSWS)中均上调表达，表达量(FC)变幅为1.23～12.24，2.89～11.37；AP2-EREBP(Unigene2016)、AP2-EREBP(Unigene442)和AP2-EREBP(Unigene160)在干旱胁迫组(LWS、 MWS和SWS)中下调，在复水组(RSWS)中又上调表达。表达量变幅为-4.65～12.22，-8.31～7.38，-6.99～6.85。bHLH(CL2529)和E2F-DP(CL4870)在干旱胁迫组(LWS、 MWS和SWS)中下调表达，变化幅度分别为-14.17～-3.73，-11.89～-4.56；在复水组中却不表达。RWP-RK(CL8167)、DBP(CL6537)，在干旱胁迫组(LWS、MWS和SWS)和复水组(RSWS)中均上调表达。BZIP(Unigene23554)在干旱胁迫组(LWS、 MWS和SWS)和复水组(RSWS)中均下调表达。综合统计数据结果显示，AP2-EREBP、NAC、E2F-DP、bHLH 和BZIP等5类转录因子对干旱和水分刺激应答显著(表 4)。

A.CK-vs-LWS的差异表达基因； B.CK-vs-MWS的差异表达基因；C.CK-vs-SWS的差异表达基因； D.RSWS-vs-SWS的差异表达基因。

表4 对比组差异表达转录因子表达量分析

注：Log2(FC)≥1时，基因表达显著。

Note:Gene expression is significant when Log2(FC)≥1.

2.4 Real-time PCR 验证

从测序数据中随机选取了9个转录因子编码基因，用Real-time PCR 方法验证测序数据的有效性。图2显示，干旱和水分刺激应答相关的AP2-EREBP、WRKY、bHLH、MYB、C2C2-Dof、RWP-RK、NAC、C2H2和ABI3VP1编码基因在处理前后表达量变化均极显著。干旱胁迫组和复水组测序数据与用荧光定量PCR 法得到的基因表达数据之间的相关性分别为0.947 8(P<0.01)，0.967 2(P<0.01)，说明干旱胁迫和复水处理后产生的差异表达基因数据有效。

均值±s，样本均为3个重复；不同大写字母表示统计结果在1%水平上有显著差异。

3 讨论与结论

本研究通过对干旱胁迫和复水处理下梭梭转录因子应答的分子机理进行研究，发掘其中的关键转录因子，对培育抗旱能力强、水分利用效率高的新作物品种具有重要的意义。研究表明，转录因子(Transcription factors，TFs)受胁迫刺激诱导产生，并将信号传递、放大，能够通过调控靶基因转录效率从多个层面降低胁迫对植物的伤害，对植物在逆境下的生长发育起到了至关重要的调节作用，尤其是AP2-EREBP、 E2F/DP、MYB、NAC、WRKY、bHLH和BZIP等家族，是信号转导通路上的关键调控因子，对干旱胁迫响应和水分刺激应答发挥了重要的调节作用[19-23]。本试验测序结果中，各对比组中 AP2-EREBP、MYB、bHLH、NAC、WRKY和ABI3VP1等转录因子编码基因的数量占主要比例，说明这几类转录因子对梭梭干旱胁迫和水分刺激应答调控作用极显著，其中，AP2-EREBP、WRKY、MYB和bHLH是所有56个家族中差异表达转录因子成员最多的蛋白家族。AP2-EREBP转录因子要调节植物对低温、干旱及高盐[23]等的分子应答反应。本研究中 AP2-EREBP是干旱和复水处理下数量丰富、表达量变化最显著的一类转录因子，且AP2-EREBP(Unigene2016)、AP2-EREBP(Unigene442)和AP2-EREBP(Unigene160)，在干旱胁迫下显著下调，而在复水组中显著上调，与 bHLH 转录因子类似，干旱胁迫处理下诱导更多的下调ERF转录因子表达，说明这类转录因子发挥调控作用的模式可能与bHLH类似。本研究结果还表明，干旱处理和复水处理后下调的bHLH编码基因数量均大于上调数，推测bHLH编码基因主要通过下调其表达调控机体应答胁迫过程。张子佳等[11]对水稻 22 个bHLH基因在PEG 和 ABA 胁迫下的表达谱进行了分析，结果表明，其中一些bHLH基因在响应干旱胁迫过程中受ABA负调控。此外，不同的bHLH基因参与应答环境胁迫的调节方式不同，表明参与 PEG 和 ABA 胁迫应答的bHLH基因具有不同的分子路径或模式。因此，梭梭中bHLH转录因子的进一步研究已在该课题组将继续进行。

Liao等[24]从研究的大豆基因中获得了 156 个 MYB 家族基因，利用酵母单杂交的方法筛选出了40个左右与逆境相关的基因，将这些基因转入拟南芥进行功能验证，发现转 GmMYB177 提高了拟南芥耐旱性。丁震乾等[25]在陆地棉中克隆了 MYB 转录因子基因GhRAX3，发现其表达在干旱胁迫 0.5 h 后即表现为显著上调，且 48 h 内持续高表达水平；而抑制GhRAX3的表达，则加快了棉花植株的失水率、细胞质膜过氧化和细胞受损程度，降低了棉花对干旱胁迫的耐受性。目前，发现的大多数 MYB 类转录因子基因均表现为提高植物抗旱能力，但亦有例外，如 Zhou 等[26]的研究任务，将JcMYB001转入拟南芥后，其过表达反而增加了拟南芥对干旱的敏感性。本研究中的MYB转录因子在干旱胁迫组与对照组相比，其下调差异表达基因较其他转录因子多，MYB家族中上调基因数随着干旱胁迫程度的增加而增加。说明本研究结果也证实了MYB家族的转录因子对干旱胁迫具有正调控和负调控性。

综合以上结果推测，发挥主要调控作用的转录因子家族或成员间可能存在一定的共表达性或相互关联性，各转录因子间协同作用对干旱和水分刺激做出精细应答调控。植物对干旱胁迫和水分刺激应答过程受众多分子和细胞通路的调控[27]。转录因子能特异与顺式作用元件结合激活或抑制下游靶基因表达，参与胁迫信号传递网络，其“阀门作用”的重要性不言而喻。本研究借助高通量测序技术，从整体层面概要性的分析了干旱和干旱后复水处理下梭梭幼苗转录因子种类及初步的表达规律，结果将为干旱胁迫应答调控网络、有目的地筛选关键调控基因和分子设计育种等研究提供基础理论数据。