干旱胁迫下谷子的转录组分析

王庆国，李臻，管延安，刘炜，潘教文

(1.山东省农业科学院生物技术研究中心/山东省作物遗传改良与生理生态重点实验室，山东 济南 250100；2.山东省农业科学院作物研究所，山东 济南 250100)

谷子[Setariaitalica(L.)Beauv.]又称为粟，是起源于我国的传统优势作物，也是粮饲兼用作物。谷子去皮后俗称小米，小米营养丰富且各种成分均衡，具有很高的营养保健价值，被我国列为“五谷”之首。谷子具有抗旱、耐瘠薄、生育期短、C4高光效、适应性广等特点，且细胞壁较厚、叶面积小、根系发达、水分利用效率较高，对干旱和盐碱胁迫等都具有较好的耐受性。随着谷子基因组测序计划的完成及基因组序列的释放，谷子逐渐发展成为研究禾本科作物新的模式作物[1,2]。

近年来结合高通量测序研究谷子的抗旱机制，并鉴定干旱响应基因及非编码RNA已有报道。Qi等[3]研究发现，PEG处理谷子后可导致2 824个基因及215个 miRNA的表达发生变化。Yadav等[4]对两个抗性不同的谷子品种测序，共鉴定到55个已知miRNA和136个新miRNA，其中18个已知和33个新miRNA均参与干旱胁迫响应。Wang等[5]在干旱敏感品种An04中鉴定到81个已知的miRNA，其中有14个表达受干旱诱导，4个表达下调，同时也鉴定到72个新的miRNA，并结合降解组鉴定到82个miRNA的靶基因。最近，Tang等[6]对两个抗性不同的谷子品种进行了生理和转录组分析，发现干旱胁迫后抗性品种 “豫谷1号”通过抑制能量代谢来保护细胞，从而维持正常的植株生长；而干旱敏感品种An04受胁迫后会消耗更多的能量，光合效率和细胞生长受到严重抑制。结合抗旱相关的QTL，作者进一步鉴定到萌发期和苗期干旱抗性相关的20个候选基因。定量蛋白质组学分析发现，谷子在干旱胁迫后会产生大量的参与渗透保护和抗氧化的蛋白，参与能量代谢的蛋白也有较明显的积累[7]。

随着谷子基因组序列的释放，谷子中胁迫相关的基因家族也逐渐被鉴定出来。如NAC、AP2/ERF、C2H2型锌指转录因子、MYB、WRKY等转录因子家族，及RNA聚合酶基因家族、PPR蛋白基因家族、DCPK基因家族、HSP基因家族等。谷子中各种胁迫及发育相关基因的功能也逐渐被研究及报道[8-14]。在水稻中异源表达自噬相关基因(autophagy-associated genes, ATG)SiATG8a能有效提高水稻对氮饥饿的耐受性[15]。在拟南芥中过表达SiCDPK24可显著增强对干旱胁迫的抗性。SiNAC1通过促进ABA的合成加速了叶片的衰老[13]。在K+缺乏和盐胁迫下，SiHAK1作为K+吸收及转运蛋白维持胞内K+的稳态。过表达LEA蛋白SiLEA14显著提高了转基因拟南芥和谷子的非生物胁迫抗性[16]。SiAGO1(argonaute1)对谷子正常的生长发育和胁迫响应是必需的，该基因缺失导致多重发育缺陷，如植株矮小、叶片狭窄卷曲、小穗以及结实率低[17]。SiLP1编码一个WRKY转录因子，该基因参与调控穗的发育、茎的伸长及种子的大小[11]。

鉴于谷子抗旱机制的研究还不够深入，本研究利用RNA-seq 技术对干旱胁迫处理不同时间的谷子进行了测序分析，筛选干旱胁迫响应基因，系统分析谷子的抗旱响应机制，为进一步克隆抗旱相关基因提供了重要依据。

1 材料与方法

1.1 谷子种植与干旱处理

精选“豫谷1号”种子播种于花盆中，每个花盆中装有3 kg营养土和沙子的混合物(比例为1∶1),在温室中正常生长2周左右，光周期为14 h，昼夜温度为30℃/25℃。用20%的PEG模拟干旱进行短期处理，将花盆中的谷子幼苗取出，用清水将根部土壤洗净，置于Hoagland营养液中适应几个小时，然后用Hoagland营养液配置20%的PEG溶液处理6 h。长期处理直接在花盆中处理，当花盆中的土壤含水量低于30%后持续处理一周的时间，使含水量控制在20%～30%。处理完的材料经液氮迅速冷冻，保存于-80℃超低温冰箱备用。

1.2 利用RNA-seq技术筛选参与谷子干旱响应的基因

与北京诺禾致源科技股份有限公司进行合作，进行了RNA-seq相关试验及数据分析，结合GO注释(包括细胞组成、分子功能、生物过程)及KEGG代谢通路注释，对筛选出的差异基因进行分析及归类。

2 结果与分析

2.1 不同时间处理后差异基因分析

在 RNA-seq 技术中，FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) 是每百万测序碱基中来自某一基因每千碱基长度的fragments数目，其同时考虑了测序深度和基因长度对fragments计数的影响，是目前最为常用的基因表达水平估算方法。利用log2FoldChange: log2(样品1的FPKM值/样品2的FPKM值)代表差异基因的变化倍数。以log2FoldChange值大于2定义为表达上调的基因，以该值小于-2为表达下调的基因。如图1A所示，与对照相比，谷子处理6 h后上调的基因有430个，下调的基因有183个。处理6 d后，与对照相比，上调的基因有606个，下调的基因有1 038个(图1B)。谷子处理6 d后与处理6 h相比，上调的基因有206个，下调的基因有611个(图1C)。

A：谷子处理6 h后的差异基因分布；B：谷子处理6 d后的差异基因分布；C：谷子处理6 d后与处理6 h相比的差异基因分布；筛选阈值默认设置为Q-value < 0.05。

对这三组差异的转录本进行共有与特有的差异转录本分析发现，有45个基因在这三组数据中共同存在，有177个差异基因在处理6 h后的谷子中被特异调控，有834个差异基因特异存在于处理6 d后的谷子中，有222个差异基因特异存在于处理6 d和6 h谷子的差异基因中(图2)。

图2 各差异转录本比较组合之间共有与特有数目的维恩图

2.2 差异基因GO (Gene Ontology)富集分析

根据筛选到的差异基因，研究特异基因在 Gene Ontology 中的分布状况将有助于阐明试验中样本差异在基因功能上的体现。在GO富集分析中，主要分为三大类biological process (BP)、cellular component (CC)和molecular function (MF)(图3)。

谷子处理6 h后的差异基因聚类分析发现，在biological process中，差异基因主要富集在biological_process和metabolic process中；在cellular component (CC)中，差异基因主要富集在extracellular region中；在molecular function (MF)中，差异基因主要富集在catalytic activity、metal ion binding和cation binding(图3A)。在干旱处理6 d后的差异基因中，在biological process中，差异基因主要富集在biological_process和metabolic process中；在cellular component (CC)中，差异基因主要富集在membrane、intrinsic to menbrane和integral to membrane；在molecular function (MF)中，差异基因主要富集在catalytic activity和oxidoreductase activity中(图3B)。在谷子处理6 d和6 h后的差异基因中，在biological process中，差异基因主要富集在metabolic process和biological_process中；在cellular component (CC)中，差异基因主要富集在extracellular region中；在molecular function (MF)中，差异基因主要富集在catalytic activity中(图3C)。

A：处理6 h；B：处理6 d；C：处理6 d和6 h。

2.3 差异基因KEGG代谢途径分析

在生物体内，不同基因相互协调行使其生物学功能，通过KEGG (kyoto encyclopedia of genes and genomes) Pathway 显著性富集能确定特异基因参与的最主要生化代谢途径和信号转导途径。通过对差异基因的KEGG途径分析发现，谷子处理6 h后的差异基因主要富集在Phenylalanine metabolism、Phenylpropanoid biosynthesis、Protein processing in endoplasmic reticulum和Biosynthesis of secondary metabolites途径中(图4A)。而处理6 d后、处理6 d和6 h的差异基因主要富集在Phenylalanine metabolism、Phenylpropanoid biosynthesis、Metabolic pathway和Biosynthesis of secondary metabolites途径中(图4B和图4C)。

A：处理6 h；B：处理6 d；C：处理6 d和6 h。纵轴表示Pathway名称，横轴表示Rich factor (候选基因集所占背景基因的比例)，点的大小表示此Pathway中差异表达基因个数多少。

3 讨论与结论

干旱会导致植物细胞脱水及渗透胁迫，严重影响植物体内各种正常的生理及代谢过程。胁迫信号会被细胞表面的受体所感知，并启动相关的信号转导途径[18-20]。在本研究中，谷子干旱胁迫不同时间后，我们鉴定到各种不同的类受体蛋白激酶(如Si030010m、Si029783m、Si000257m)的表达明显上调。干旱胁迫会导致胞内Ca2+浓度增加，钙信号被相关的Ca2+结合蛋白(如Si002992m、Si024895m、Si002992m)感知后启动钙信号转导。胁迫信号通过激活一些蛋白激酶或者磷酸酶进一步向下传递[18]。促分裂原活化蛋白激酶(mitogen-activated protein kinase, MAPK)级联途径是真核生物中普遍存在的信号转导模式，渗透胁迫能够激活MPK3、MPK6和MPK4，干旱胁迫能够激活AtMEKK1-MKK1-MPK4级联途径[21,22]。在谷子中我们同样鉴定到一个MAPKKK(Si025330m)基因被明显诱导激活。对于特定刺激的细胞响应主要依赖于MAPK激活的响应幅度和持续时间，因此这个过程受到精细的调控。蛋白磷酸酶是MAPK信号转导途径重要的负调控因子。PP2C型的磷酸酶AP2C1能够灭活MPK4和MPK6[23]。在谷子干旱胁迫后有10个PP2C型蛋白磷酸酶的表达被激活，从而精细调控干旱胁迫信号转导。

激活后的MAPK通过磷酸化下游的转录因子进而调控胁迫基因的表达。胁迫激活的MPK3/MPK6/MPK4能够磷酸化WRKY和MYB转录因子[21]。最近的报道也发现，类受体蛋白激酶GUDK(GROWTH UNDER DROUGHT KINASE)磷酸化并激活转录因子OsAP37介导干旱胁迫信号转导，激活后的OsAP37进一步诱导胁迫响应基因的表达来增强抗旱性和水稻产量[24]。本研究谷子在干旱胁迫后，多种转录因子(如DREB、NAC、WRKY、MYB、bHLH等)的表达被明显诱导，这些结果与近两年来谷子中转录因子家族鉴定的结果基本一致。激活后的转录因子通过结合胁迫响应基因启动子中的顺式调控元件，如DRE/CRT、ABRE、YCRS/MYBRS，进而调控基因的表达[18]。下游胁迫响应的功能蛋白，如LEA(late embryogenesis abundant protein)、不同分子量的HSPs(heat shock protein)、脱水素蛋白dehydrin、分子伴侣蛋白ClpB1(chaperone protein)等在干旱处理6 h和6 d后被明显诱导。进一步分析发现，水孔蛋白(aquaporin)在干旱胁迫6 d后明显下调，水孔蛋白的下调可能有助于保持植物体内的水分，防止水分的散失。泛素蛋白酶体途径是真核生物体内重要的蛋白质降解途径，E3泛素连接酶在该过程中能特异性识别并结合底物蛋白，从而调节靶蛋白的降解代谢。研究表明E3泛素连接酶在植物胁迫响应过程中起着重要的调控作用。拟南芥中的E3泛素连接酶DRIP1和DRIP2通过结合并泛素化DREB2A抑制其调控的干旱胁迫响应基因的表达[25]。在本研究谷子胁迫响应的后期(干旱处理6 d后)，发现10个E3泛素连接酶基因的表达受到明显调控，说明E3泛素连接酶介导的泛素化途径在谷子胁迫响应过程中起着重要的调控作用。

本研究对干旱胁迫下的谷子进行了转录组分析，比较了不同时间处理后的差异表达基因，为进一步筛选并克隆谷子关键抗旱基因提供了依据，也为揭示谷子的抗旱机制奠定了基础。