基于RNA-seq的谷子萌芽期抗旱相关基因挖掘与分析

代小冬 朱灿灿 宋迎辉 王春义 代书桃 秦 娜 李君霞

(河南省农业科学院粮食作物研究所，河南郑州 450002)

干旱是限制作物生长发育的主要非生物逆境灾害[1]。随着气候变化和世界人口的不断增长，水资源短缺的问题更加突显，干旱发生的频率和严重程度将会进一步加剧[2]。因此，解析作物抗旱机制、培育抗旱作物品种对减轻干旱灾害及其造成的损失具有重要的意义。

谷子[Setaria italica(L.) Beauv]起源于中国，迄今已有8 700 多年的栽培历史[3]。谷子是二倍体自花授粉植物，且拥有较小的基因组(约430 Mb)。谷子具有生育期短、单株籽粒多、适应性广、抗旱、耐瘠薄等特点，使其渐渐成为抗旱研究的模式作物。在抗旱性鉴定方法方面，前人相继确立了谷子萌芽期[4－5]、孕穗期[6]和全生育期[7]的抗旱性鉴定指标，并进行了抗旱性鉴定。在抗旱相关QTL 定位方面，Qie 等[8]利用豫谷1 号与青狗尾草构建的重组自交系(recombinant inbred lines，RIL)群体定位到18 个与谷子抗旱相关的QTL。在参与谷子抗旱的功能基因与蛋白挖掘方面，Tang 等[9]对抗旱品种豫谷1 号和干旱敏感品种安04 在干旱条件下的转录组进行分析，并结合抗旱相关QTL 定位结果，确定了20 个抗旱候选基因。赵晋锋等[10]研究认为SiCIPK19 基因参与了谷子拔节期、抽穗期和灌浆期对干旱胁迫的响应，推测该基因参与谷子对非生物逆境的应答，尤其在抽穗期和灌浆期干旱胁迫应答中发挥重要作用。Wang 等[11]通过对谷子中149 个bHLH 转录因子的研究，发现其可能通过调控下游抗旱基因的表达参与谷子的抗旱反应。Feng 等[12]通过分析谷子中39 个核转录因子的功能，发现SiNF-YA1 和SiNF-YB8 可以激活相关抗旱基因，提升谷子生理状态，提高谷子抗旱能力。Pan 等[13]利用定量蛋白质组学分析的方法发现了321个与谷子抗旱性相关的蛋白。

目前虽然在谷子抗旱性研究方面取得了一定的进展，但是由于谷子抗旱机制复杂和研究起步较晚，对谷子抗旱分子调节机制的认识还比较有限，因此在谷子抗旱关键基因发掘和干旱分子响应机制方面的研究有待进一步加强。本研究以谷子抗旱品种山西2010 和干旱敏感品种K359*M4－1 为材料，应用转录组测序RNA-seq 技术对两个品种干旱胁迫前后萌发期种子进行转录组测定，获得抗旱相关差异表达基因(differentially expressed genes，DEG)，并通过GO(gene ontology)富集和KEGG(kyoto encyclopedia of genes and genomes)分析对DEG 进行功能分类，以期挖掘谷子抗旱相关基因，为解析谷子抗旱遗传机制和指导谷子抗旱育种提供理论依据。

1 材料与方法

1.1 试验材料

供试材料K359*M4－1 和山西2010 分别来源于河北省农林科学院谷子研究所和山西省农业科学院经济作物研究所。

1.2 谷子萌芽期干旱胁迫处理

选用0、10%、15%、20%、25%、30% 6 个聚乙二醇6000(polyethylene glycol 6000，PEG－6000)浓度进行预试验，筛选出谷子萌芽期抗旱鉴定的PEG 6000 适宜浓度为20%。利用浓度为20%的PEG 6000 对山西2010 和K359*M4－1 进行萌芽期抗旱性鉴定。每份材料选择饱满的种子，首先用5%次氯酸钠消毒15 min，然后用ddH2O 冲洗干净。选择50 粒灭过菌的种子放入直径为9 cm 的培养皿中，再将培养皿放入培养箱，28℃暗培养。干旱胁迫组培养皿中加入5 mL 20%PEG 6000，对照组培养皿中加入5 mL ddH2O。每个试验设3 次重复。

分别调查第2、第4、第6、第8 天种子的发芽数，并将第2、第4、第6、第8 天种子的萌发率记为nd2、nd4、nd6、nd8。参考Bouslama 等[14]的方法计算萌发指数(promptness index，PI)和萌发抗旱指数(sprout index of drought resisting，SIDR)，PI ＝1.00nd2 ＋0.75nd4 ＋0.50nd6 ＋0.25nd8，SIDR＝处理萌发指数/对照萌发指数。

1.3 RNA 文库构建和测序

文库构建和转录组测序由北京诺禾致源生物信息科技有限公司完成。

1.4 数据的处理与筛选

序列数据由Illumina HiSeq 2500 测序得到的原始图像数据经base calling 转化后形成FASTQ 格式的文件，该序列数据包含reads 的序列及碱基的测序质量，在经过去除和过滤含adapter、N 过多或低质量碱基的dirty raw reads 后获得高质量的clean reads，用于后续的信息分析。

利用转录组数据比对软件TopHat2 将所有的clean reads 与谷子(豫谷1 号)参考基因组序列进行比对。利用比对的reads 计算基因的表达量以确保结果的高准确性。用FPKM (fragments per kilobase of transcript per million mapped reads)值衡量计算表达量。分别利用GOseq R 和KOBAS 软件包进行差异表达基因的GO 和KEGG 富集分析(伪发现率false discovery rate，FDR<0.05)。

本研究报道的原始数据已存入中国科学院北京基因组研究所大数据中心基因组序列档案库，accession ID:CRA002537。公众可访问网址:https:/ /bigd.big.ac.cn/gsa 获取数据。

1.5 实时荧光定量PCR 分析

为了验证谷子萌芽期干旱胁迫转录组结果的准确性，从DEG 数据库中随机挑选28 个基因，利用Primer 6.0 设计引物(表1)，分析实时荧光定量PCR(realtime quantitative PCR，RT-qPCR)结果是否与转录组测序结果一致。

表1 RT-qPCR 引物序列Table 1 RT-qPCR primer sequence

2 结果与分析

2.1 干旱处理对谷子发芽的影响

以浓度为20%的PEG 6000 对91 份谷子材料的萌芽期抗旱性进行评价。由图1可知，山西2010 的萌发抗旱指数(sprout index of drought resisting，SIDR)达到92%，抗旱性较强，而K359*M4－1 的萌发抗旱指数仅为64%，抗旱性较弱。

根管充填即刻质量评价标准[1]（1）适充：根充材料距解剖根尖≤2mm,根管三维充填封闭严密；（2）超充：根管充填材料超出根尖孔；（3）欠充：根充材料距解剖根尖＞2mm或者根管三维充填不严密，X线显示根充物不致密。

图1 干旱胁迫对山西2010 和K359*M4－1 发芽的影响Fig.1 Effects of drought stress on germination of Shanxi 2010 and K359*M4－1

2.2 干旱处理萌芽期谷子种子的转录组分析

利用RNA-seq 技术对山西2010 和K359*M4－1萌芽期种子进行转录组测序。将干旱处理前后山西2010 样品分别命名为SC 和ST，将干旱处理前后K359*M4－1 样品分别命名为KC 和KT。将校正后P<0.05 作为标准筛选每个比较组的DEG。结果显示，山西2010 干旱处理前后(SC vs ST)共鉴定到2 300 个DEG，其中上调基因有1 002 个，下调基因有1 298 个。K359*M4－1 干旱处理前后(KC vs KT)共鉴定到3 652个DEG，其中上调基因1 861 个，下调基因1 791个。综上可知，在干旱胁迫下，与抗旱品种山西2010相比，干旱敏感品种K359*M4－1 的代谢发生了更大的变化。

韦恩图分析显示，干旱胁迫处理前后山西2010 和K359*M4－1 之间存在1 268 个共有DEG，在山西2010 和K359*M4－1 2 个品种中分别存在1 032 个和2 384 个特有的DEG(图2)。这些基因参与了谷子对干旱胁迫的响应，并可以用来解释2 个谷子品种的抗旱性差异。

图2 DEG 表达韦恩图Fig.2 Venn diagrams of DEGs among different comparisons

2.3 DEG 功能分析

为了进一步探索谷子对干旱胁迫的响应机制，将矫正后P<0.05、阈值(｜Log 2 fold change ｜≥2)作为标准筛选每个比较组的DEG。分别在山西2010 和K359*M4－1 中鉴定到82 个(40 个上调基因，42 个下基因调)和112 个(51 个上调基因，61 个下调基因)DEG(表2)。DEG 主要为编码类锌诱导的促进因子蛋白、类萌发素蛋白、蛋白质磷酸酶、转运蛋白、胚胎发育晚期丰富蛋白、转录因子、过氧化物酶等的基因。

表2 部分谷子萌芽期响应干旱胁迫DEGTable 2 Part of DEGs response to drought stress at germination stage of foxtail millet under drought stress

2.4 DEG 功能富集分析

通过对不同比较组鉴定到的DEG 进行GO 和KEGG 代谢途径富集分析。根据GO 功能注释，SC vs ST 和KC vs KT 中的DEG 分别富集在52 个和21 个生物学过程。主要集中在催化活性(GO:0003824)、单个有机体代谢过程(GO:0044710)、氧化还原过程(GO:0055114)、血红素结合(GO:0020037)、四吡咯结合(GO:0046906)等生物学过程(表3)。

表3 干旱胁迫下谷子DEG 的GO 富集分析Table 3 GO enrichment analysis of s DEGs in foxtail millet under drought stress

KEGG 富集分析表明，不同比较组鉴定到的DEG主要富集在淀粉和蔗糖代谢，植物激素信号转导，光合作用－天线蛋白，苯丙素生物合成，角质、亚氨酸和蜡生物合成，次生代谢物的生物合成等途径(图3)。

图3 干旱胁迫下山西2010(A)和K359*M4－1(B)中DEG 的KEGG 富集图Fig.3 KEGG enrichments map of the DEGs in Shanxi 2010(A) and K359*M4－1 (B) under drought stress

2.5 实时荧光定量PCR 验证

图4 实时荧光定量PCR 结果Fig.4 Confirmation of the transcriptomic profiles of selected genes by RT-qPCR

3 讨论

RNA-seq 可以全面地获得某一物种特定细胞或组织在某一状态下的几乎所有的转录本及基因序列，是一种快速、有效鉴定差异表达基因的方法，是揭示植物抗逆机理以及筛选抗逆基因的主要研究技术[15－20]。本研究利用RNA-seq 在全基因组表达水平上鉴定了谷子萌芽期的抗旱基因。通过对干旱处理前后抗旱品种和干旱敏感品种的转录组比较，在山西2010 和K359*M4－1 中分别发现了2 300 个和3 652 个可能参与干旱胁迫响应的DEG。这些基因的表达和调控能使谷子抵御干旱胁迫的危害。进一步比较发现，K359*M4－1 中的DEG 明显多于山西2010 中的DEG，表明在干旱胁迫下，与抗旱品种相比，干旱敏感品种的代谢发生了更大的变化。上述结果与前人研究结果一致[21－22]。

本研究中，DEG 功能富集分析结果表明，DEG 显著富集于氧化还原相关的生物学过程，表明植物可以通过调控活性氧代谢来保护细胞膜的稳定性，从而提高植物的抗旱性[23－24]。苯丙素生物合成、次生代谢物生物合成、淀粉和蔗糖代谢、苯丙氨酸代谢等途径显著富集在两个比较组中。植物次生代谢产物是在长期进化过程中植物与生物和非生物因素相互作用的结果，在植物对环境胁迫的适应过程中发挥重要作用[25]。当植物处于干旱环境时，植物次生代谢产物会通过减缓淀粉代谢来保存能量以延长细胞的存活时间[15]。说明富集在这些途径上的基因在谷子响应干旱胁迫方面发挥了重要作用。

本研究获得了很多与干旱胁迫相关的DEG，如编码类锌诱导促进因子蛋白的基因Seita.3G089200 在两个比较组中的表达都表现明显上调。前人研究表明，类锌诱导促进因子可以通过调节气孔的关闭提高植物的抗旱性[26]，因此，推测Seita.3G089200 参与了谷子对干旱的响应。

在植物中，钾不仅是生长发育所必需的宏量营养素，而且是介导植物响应胁迫的主要信号[27]。钾吸收透性酶/高亲和力钾转运体/钾转运体(K＋uptake permease/high-affinity K＋transporter/K＋transporter，KUP/HAK/KT)转运蛋白在调控甘蔗响应低钾、干旱和盐胁迫方面发挥了重要作用[28]。本研究发现2 个编码KUP/HAK/KT 转运蛋白的基因(Seita.9G182400和Seita.9G182500)在干旱胁迫下上调表达。表明这些基因可能是参与谷子抗旱的关键基因。ABC 转运蛋白参与跨膜运输、信号传导和细胞解毒等多种重要的生理过程，使植物能够适应环境的变化，并调节植物对生物和非生物胁迫的反应[29]。在拟南芥(Arabidopsis thaliana)中超表达AtABCG36 可以显著增强相株的抗旱性[30]。

本研究发现4 个编码ABC 转运蛋白的基因(Seita.9G001100、Seita.9G447200、Seita.3G042600 和Seita.2G290000)参与了谷子对干旱胁迫的响应。类萌发素蛋白基因受多种环境信号诱导，在植物生长发育及对非生物胁迫、生物胁迫的响应中起着重要的作用[31－34]。超表达大豆(Glycine max)的GmGLP7 可以提高拟南芥的抗旱性和抗氧化能力[35]。本研究有10个编码类萌发素蛋白的基因在干旱胁迫下差异表达，推测这些基因与谷子抗旱性相关。蛋白质磷酸酶是许多信号转导途径的调控因子，在植物抗氧化、耐冷、抗旱等过程中起关键作用[36]。如烟草(Nicotiana tabacum)中的NtPP2C1 的表达受干旱诱导，在干旱胁迫条件下的表达显著提高[37]。本研究中，在干旱条件下，谷子中3 个编码蛋白质磷酸酶的基因(Seita.7G021400、Seita.6G005300 和Seita.3G218800)表达上调，可能参与了谷子对干旱胁迫的应答。转录因子可通过调控多个下游靶基因的表达，诱导多种非生物胁迫的保护机制[38－40]。如NAC、MYB、bHLH、AP2/ERF 等转录因子直接或间接参与抗旱响应[41－43]。本研究发现编码NAC 转录因子基因Seita.2G028200、编码MYB 相关蛋白基因Seita.1G374700、bH035 的同源基因Seita.7G045200 的表达在干旱胁迫前后差异显著，推测这些基因可能通过调节相关基因表达参与谷子对干旱的响应。LEA 蛋白参与植物对各种非生物和生物胁迫响应，包括盐害[44]、干旱[45]、热害[46]、冷害[47]等。提高LEA 蛋白基因的表达可以提高植物的耐旱性[48－50]。本研究中，编码LEA 蛋白的Seita.2G437900、Seita.9G420500、Seita.5G287800 和Seita.3G162400 在干旱胁迫下表达显著上调，表明这些基因可能与谷子抗旱相关。当植物遭受干旱胁迫时，其通过抗氧化酶系统保护组织免受由于活性氧积累导致的氧化损伤[51]。过氧化物酶是一种抗氧化酶，在植物抗旱胁迫过程中发挥重要的作用[52]。本研究发现，Seita.7G247400 和Seita.3G052900 在谷子受到干旱胁迫时，表达显著上升，推测它们可能参与了干旱胁迫下谷子中活性氧的清除过程。

4 结论

本研究通过对干旱胁迫前后2 个谷子品种的萌芽期种子进行转录组分析，获得了一批包括编码类萌发素蛋白、蛋白磷酸酶、LEA 蛋白、转录因子、ABC 转运蛋白等与干旱胁迫相关的基因。DEG 显著富集在淀粉和蔗糖代谢，植物激素信号转导，光合作用－天线蛋白，苯丙素生物合成，角质、亚氨酸和蜡生物合成，次生代谢物生物合成等途径，表明这些途径参与了谷子的抗旱反应。本研究结果为解析谷子抗旱机制、筛选谷子抗旱品种奠定了很好的理论基础。