PEG胁迫下玉米转录组差异性分析

124

(1.河南牧业经济学院食品与生物工程学院(酒业学院),郑州 450046;2.河南牧业经济学院河南省白酒风格工程技术研究中心,郑州 450046;3.河南牧业经济学院郑州市白酒酿造微生物技术重点实验室,郑州 450046;4.张弓老酒酒业有限公司,河南 宁陵 476733)

玉米是世界上最重要的粮食作物之一，不仅可食用，而且还广泛用于白酒酿造及酒精生产；但其又是旱地作物中对水需求量较多，并且对水分胁迫较为敏感的作物[1-3]。在中国，玉米有着悠久的种植和培育历史，其种植面积约占耕地面积的25%，年产量可达1.2亿t。而干旱是限制玉米生产发展的第一要素，导致70%以上的玉米遭受干旱威胁，每年造成的产量损失总计达到1 500万t以上[4-7]，因此研究玉米的抗旱机理尤为重要。

刘成等研究发现，在干旱条件下，各种酶活性增加并促进体内合成代谢，以增加玉米抗旱性[8]。贾晓燕等基于生理生化指标研究得出，细胞中SOD、POD、CAT、MDA、脯氨酸积累和丙二醛含量变化与作物抗旱性相关的结论[9]。Ali Noman等发现，干旱胁迫增强了玉米植株SOD和POD的活性以改善玉米植株的抗氧化防御系统，减轻干旱胁迫对玉米植株的危害[10]。前人对玉米抗旱性研究多涉及生理生化指标，但较少涉及从RNA水平研究玉米抗旱基因的表达情况。

而随着后基因组时代的到来，转录组测序技术(RNA-seq)广泛应用于医学[11]、药物研发[12]和农业科学[13]等基础研究[14]。对玉米转录组序列的结构进行解释、对转录组的功能进行分析，能从分子水平揭示玉米植株内部的生物抗性[15]。目前RNA-seq应用于玉米相关性研究较多；谷荣欢利用RNA-seq发现玉米K 305 ms花药与发育相关的代谢网络中部分基因表达量不足，是导致雄配子形态和发育过程异常的原因[16]。李川等报道，花粒期的玉米在高温胁迫下，热激蛋白、热激转录因子、生长素基因和磷脂酰基醇-4,5-二磷酸基因显著表达[17]。王其等研究发现，将杜仲内生拮抗细菌DZSY 21引入玉米，可调节玉米叶片中抗性相关基因的表达，适应生存环境的变化，减少伤害[18]。

但是利用转录组测序分析对比PEG胁迫前后玉米基因差异表达的研究较少报道。因此本研究在RNA-seq的基础上，通过对玉米转录组进行一系列的数据比对、分析，初步解释玉米抗旱应答机制，以期为玉米的抗旱性研究及玉米抗旱基因的发掘奠定了基础。

1 材料与方法

1.1 实验材料

试验材料为玉米自交系郑58。

1.2 实验设计

采用人工控水盆栽试验法：将玉米移栽在花盆中(花盆直径30 cm，高25 cm，每盆装15 kg培养基质)生长，盆中含有培养基质(3园土∶1珍珠岩∶1泥炭土)，一次性施入3 g复合肥(N∶P2O5∶K2O=15∶9∶10)，种前保持土壤含水量不低于85%以保证出苗。设置样品S 1处于正常供水和S 3处于干旱胁迫2个状态，3次重复，对照组玉米保持整个生长期正常浇水的湿润状态，而干旱处理从苗期(5叶期)开始控制灌溉20%的PEG-6000渗透调节剂模拟干旱胁迫[19]。

1.3 RNA提取

改进的SDS酚法取：

1) 取0.1 g新鲜玉米叶片，液氮研磨成粉末状；

2) 转入1.5 mL离心管中，迅速加入600μL SDS提取液(2% SDS，25 mmol·L-1EDTA，100 mmol·L-1Tris-HCl，2%PVP)和40μL巯基乙醇，混匀后再加入150μL无水乙醇和150μL 3 mol·L-1的KAc，混匀；

3) 在转速12 000 r·min-1、4 ℃下离心15 min，取上清，加入250μL水饱和酚和250μL氯仿，混匀后冰上放置10 min；

4) 在转速12 000 r·min-1、4 ℃下离心15 min，取上清，加入等体积的氯仿抽提，混匀后冰上放置10 min；

5) 在转速12 000 r·min-1、4 ℃下离心15 min；取上清，加入等体积异丙醇，混匀后冰上放置10 min；

6) 在转速12 000 r·min-1、4 ℃下离心20 min，弃上清，加入75%的乙醇1 mL，温和洗涤沉淀；

7) 在转速12 000 r·min-1、4 ℃下离心5 min，弃上清，干燥后用50μL DEPC水溶解RNA[20]。

1.4 文库构建及RNA-Seq测序

样品检测合格后，先采用Oligo(dT)方法分离并纯化出mRNA，并将mRNA打断为片段，通过RT-PCR扩增得到玉米转录组cDNA文库，再采用Illumina PE 150测序平台技术对构建好的文库进行高通量测序。

表1 转录组测序数据及质量情况统计表

1.5 数据的处理

使用软件：fastx-toolkit-0.0.14对每一个样本的原始测序数据进行测序相关质量评估；使用StringTie-1.3.3 b软件对每个样本进行拼接、组装；使用软件DEGseq 2-1.24.0进行组间差异表达分析，默认参数：p-adjust<0.001 &|log 2 FC| ≥1；使用BLAST 2 GO-2.5与GO数据库进行序列比对；使用DIAMOND-0.9.24将转录本与NCBI-NR，Swiss-Prot，COG数据库进行序列比对；采用HMMER-3.2.1与Pfam数据库比对；采用软件Goatools-0.6.5进行GO功能富集分析；使用KOBAS-2.1.1获得到KEGG Orthology结果。

2 结果与分析

2.1 玉米转录组测序数据及质量情况

运用统计学的方法对所有测序reads的每个cycle进行碱基分布和质量波动的统计；从表1可知，对测序数据过滤后S 1、S 3分别得到45 426 908、48 371 948条clean reads，且每个样本的clean bases均在6.78 Gb以上；碱基百分比Q 20均大于98.26%、碱基百分比Q 30均大于94.65%，其GC含量在55.9%～56.09%之间，含量接近；碱基测序错误率在0.024 1%～0.024 5%的较低水平；以上数据直观地反映出测序样本的测序质量和文库构建质量都比较好，能用于后续的分析。

2.2 转录本组装及功能注释情况

使用StringTie软件对2个样本进行拼接，最终merge在一起，结果见图1，在该转录本长度范围内共得到161 823个转录产物(transcipts)，转录本长度主要集中在600～3 000 bp之间；利用生物信息学方法将转录本组装获得的转录本与六大数据库NCBI-NR、Swiss-Prot、Pfam、COG、GO和KEGG进行比对，全面获得转录本的功能信息并对各数据库注释情况进行统计。从表2可知，共有138 270个转录本得到注释，在e-value≤1 e-5条件下，Nr注释成功的transcript数目有133 087条，占组装得到总transcript的96.25%；GO注释成功的transcript数目有108 300条，占78.33%；Swiss-port注释成功的transcript数目有109 166条，占78.95%；COG注释成功的transcript数目有130 072条，占94.07%；KEGG注释成功的transcript数目有65987条，占47.72%；Pfaw注释成功的transcript数目有99 546条，占71.99%。

图1 转录本组装长度分布统计

表2 转录本功能注释情况统计

2.3 玉米转录组响应干旱的差异表达分析

利用DEGseq软件，参数设置：p-adjust<0.001 &|log 2 FC|≥1，以未进行干旱胁迫的玉米转录组S 1为参考，与转录组S 3进行数据比对；从图2可知，玉米转录组共获得12 358个差异表达基因，其中上调基因5 275个，下调基因7 083个。

2.4 差异表达基因的COG功能分析

对差异表达基因进行COG功能分类，结果见表3，共得到23个不同的COG功能，2 127个转录本参与细胞过程和信号传导的9个COG功能，其中翻译后修饰，蛋白质转换，伴侣蛋白681个，占7.74%；信号传导机制功能有转录本675个，占7.67%；非细胞运输、分泌和囊泡运输391个，占4.44%。1 337个转录本归纳到新陈代谢类的8个COG功能，其中碳水化合物运输和代谢315个，占3.58%；氨基酸运输和代谢265个，占3.01%。1 296个转录本参与信息存储与处理的5个功能，转录506个，占5.75%；翻译，核糖体结构和生物发生356个，占4.05%。而分类到POORLY CHARACTERIZED未知功能的转录本有4 039个(45.91%)。

表3 差异表达基因COG功能分类

图2 基因差异表达火山图分析

2.5 干旱胁迫下玉米差异表达基因的分析

2.5.1玉米转录组差异表达基因GO功能富集分析

采用软件Goatools对干旱胁迫处理后玉米的差异表达基因集中的转录本进行GO功能富集分析。结果见图3，可以看出，玉米转录组差异表达基因富集在生物过程、细胞成分、分子功能三大类；差异表达基因富集的前20亚类中生物过程大类中有6个亚类，分别是光合作用、mRNA加工、tRNA代谢过程、ncRNA代谢过程、细胞表面受体信号通路、代谢过程。还富集在细胞成分大类中的8个亚类，分别为叶绿体、质体、细胞质、质体部分、叶绿体的部分、胞质部分、细胞内的部分、细胞部分。分子功能大类的6个亚类得到富集，分别是GTP活性、核糖核酸结合、三磷酸鸟苷结合和鸟苷核糖核苷酸结合、鸟苷核苷酸结合、嘌呤核苷绑定。

图3 差异表达基因GO功能富集分析

2.5.2玉米转录组差异表达基因KEGG功能富集分析

Pvalue corrected小于0.05时认为显著富集，基于此，对玉米转录组参与KEGG代谢通路进行富集，结果见图3，玉米转录组参与KEGG代谢通路分为三大分支：代谢、遗传信息处理和细胞过程，共获得28个标准KEGG代谢通路。其中代谢分支中富集到的通路有光合作用-天线蛋白、光合生物的固碳作用、卟啉和叶绿素代谢、脂肪酸降解、精氨酸和脯氨酸代谢、磷酸肌醇信号、糖酵解和糖异生、丙氨酸代谢、氨基糖和核苷酸糖代谢、淀粉和蔗糖代谢、氰氨基酸代谢、光合作用、类固醇生物合成、牛磺酸和次牛磺酸代谢、精氨酸生物合成、花生四烯酸代谢、缬氨酸、亮氨酸和异亮氨酸的降解、硫代谢、亚麻酸代谢、缬氨酸、亮氨酸和异亮氨酸生物合成、亚油酸的新陈代谢、烟酸和烟酰胺代谢、苯丙合成。遗传信息处理富集到的通路有氨酰生物合成、蛋白酶体、蛋白输出、剪接体。细胞过程富集到的通路有吞噬体。

图4 差异表达基因KEGG富集分析

3 结果与讨论

近年来，转录组测序技术的应用很大程度上推动了植物分子水平方面的研究[21]。本研究对正常灌溉和PEG模拟干旱胁迫下的玉米进行RNA-sep测序，玉米转录组共获得12 358个差异表达基因，其中上调基因5 275个，下调基因7 083个；可见，抗旱基因为上调表达，干旱胁迫限制了玉米的某些代谢活动；下调表达高于上调表达；这与王伟东[22]、吉福桑[23]和杨文军等[24]对植物干旱胁迫下的基因表达上下调结果一致。对差异表达基因进行COG功能分类，共得到23个不同的COG功能，细胞过程和信号传导中的9个COG功能、新陈代谢类的8个COG功能、信息存储与处理的5个功能以及POORLY CHARACTERIZED的未知功能；除未知功能的转录本4 039个外；翻译后修饰，蛋白质转换，伴侣蛋白681个，占7.74%；信号传导机制功能有转录本675个，占7.67%；转录506个，占5.75%；非细胞运输、分泌和囊泡运输391个，占4.44%。碳水化合物运输和代谢315个，占3.58%；翻译，核糖体结构和生物发生356个，占4.05%；氨基酸运输和代谢265个，占3.01%。植物可通过调整碳水化合物的合成利用来满足植物对能量的需求，以缓解干旱胁迫[25]，同时干旱胁迫会增强玉米植株SOD等酶的活性以改善玉米植株的抗氧化防御系统[10]，而SOD的合成表达又由其伴侣蛋白来实现[26]。由此可见，玉米在PEG胁迫下通过信号传导环境变化，进而调控转录翻译以合成SOD等酶，提高玉米适应干旱环境的抗性。

差异表达基因GO功能富集分析，说明玉米转录组特性与生物过程、细胞成分、分子功能相关；进一步对KEGG代谢通路富集，在Pvalue corrected小于0.05时显著富集到28条标准KEGG代谢通路，这些基因主要参与光合作用-天线蛋白、光合生物的固碳作用、卟啉和叶绿素代谢、脂肪酸降解、精氨酸和脯氨酸代谢、磷酸肌醇信号、糖酵解和糖异生等生命代谢活动。其中多条代谢通路涉及光合作用等相关途径，说明干旱胁迫对光合作用产生了显著影响，刘素军等报道，干旱胁迫会加强光合作用-天线蛋白途径向相连的反应中心输送能量，从而响应水分胁迫环境[27]。孙琪研究发现，光合生物的固碳作用用于合成更多的有机物满足植物对能量的需求，以此来提高抗盐能力[28]。张国伟研究发现，脂肪酸含量是影响光合作用的重要因素之一[29]。Trovato M等发现，脯氨酸的积累是生物有机体对环境胁迫产生的一种应答反应，能够增强植物的抗逆能力[30]。王建等报道，精氨酸代谢会间接促进PA、NO的合成，从而提高番茄对Cu胁迫的抵抗能力[31]。Ishitani M等发现，磷酸肌醇信号途径在高等植物的生长、发育及对环境胁迫反应中发挥重要作用[32]。上述代谢通路都与植物的抗逆性相关，但个别代谢通路在玉米中的调控机理暂不清楚，这也为进一步了解玉米抗旱机理指明了研究方向。

RNA-Seq是一种在转录水平上更为精确的测定分析方法，其通量高、成本低、分辨率高、灵敏度高且不受物种限制[33-35]。本实验采用RNA-Seq研究玉米转录组，并进行数据分析工作，其结果极大丰富了玉米转录组数据资源，为探讨玉米抗旱机理打下了坚实的基础。