转录组及代谢组联合解析玉米响应花粒期高温胁迫机制

李 川，乔江方，黄 璐，张美微，张盼盼，牛 军，刘京宝

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

随着大气温室效应加剧，近年来我国短期异常高温频发[1]。黄淮海平原是我国夏玉米主要产区[2]。7-8月为此地区年高温期，恰逢夏玉米对温度最敏感的抽雄吐丝期。日平均温度25～28 ℃，空气相对湿度70%最有利于玉米开花授粉[3]。温度超过35 ℃不利于玉米授粉，超过39 ℃则引起花粉失活，造成热害导致减产甚至绝收[4]。高温除影响玉米雄穗的发育，使雄穗分枝变小、数量减少、小花退化、花药瘦瘪、花粉活力降低外[5]，还影响雌穗的发育。主要表现在分化异常、延缓吐丝，造成雌雄不协调，授粉结实不良[6]。高温还能引起气孔关闭[7]、破坏叶绿体结构、降低蛋白酶活性，从而减弱光合作用[8]。高温加速生育过程中各种生理反应[9]，缩短生育期。灌浆期高温缩短灌浆时间[10]，减少干物质积累量[11]、籽粒容重，最终降低产量和品质[12]。受害程度随着高温胁迫持续时间及强度增加而增强[13]，不同基因型玉米响应高温胁迫的机制不尽相同[14]。

玉米是我国重要的粮食、饲料、工业三元作物。全球气候变暖高温天气频发[15]，尤其是抽雄吐丝期高温胁迫，改变了玉米的生育进程。李树岩等[16]报道了随着高温时间增加，河南省夏玉米营养生长期平均缩短3.1～4.7 d，全生育期将缩短8.6～12.9 d。目前对玉米高温胁迫的研究多从形态、生理生化指标以及产量、品质性状开展。高温胁迫影响叶片结构及玉米形态构成。叶绿体是植物进行光合作用的重要细胞器。随着高温胁迫时间的延长，叶绿体膨大变形，叶绿体被膜断裂甚至解体，类囊体片层松散、断裂，叶绿体基质外溢，从而抑制光合作用[17-18]。营养生长是植物生长发育的基础。玉米苗期高温胁迫下叶片发育变慢、发黄甚至早衰[19]，叶面积减少，植株瘦弱[20]。抽雄期高温胁迫下根茎干物质量降低，株高降低，叶片卷曲。高温胁迫还影响抗氧化生理特征。高温胁迫条件破坏植株生长发育中活性氧产生与清除的动态平衡[21]。活性氧大量积累损伤膜系统，电解质流出，细胞生理生化特性发生紊乱。高温逆境下超氧化物歧化酶(SOD)、过氧化物酶(POD)、过氧化氢酶(CAT)等抗氧化酶及脯氨酸、可溶性蛋白、可溶性糖等渗透调节物质共同协作清除积累的过多活性氧，维持细胞功能[22]。王永军等[23]发现，玉米叶片衰老与SOD、POD活性呈负相关，不同部位叶片抗氧化酶活性不同。渗透调节物质可以保持细胞含水量，维持细胞膨压，防止细胞失水，维持其正常功能。胡秀丽等[24]把SOD作为作物耐高温的评价指标之一。渗透调节具有局限性，随着高温胁迫延长，渗透调节能力降低。高温胁迫降低作物光合作用。高温胁迫初期，为了防止水分过度蒸发，气孔关闭，胞间CO2浓度降低，光合反应受到抑制；随着高温持续，尽管为了防止高温灼伤气孔打开，但是光合蛋白酶活性降低，叶绿体结构破坏，光合作用被抑制[25]。另有研究报道，高温胁迫还可以通过降低光合系统Ⅱ(PS Ⅱ)原初电子交换效率、荧光猝灭系数、碳代谢关键酶活性而降低作物光合能力[26]。不同作物不同基因型不同生长发育阶段对高温响应机理不同。张吉旺等[8]报道，玉米生殖生长阶段高温胁迫时核酮糖-1，5-二磷酸羧化酶(RuBPCase)、磷酸烯醇式丙酮酸羧化酶(PEPCase)活性降低，气孔导度及光合速率均降低。高温胁迫影响作物内源激素。植物内源激素是一类在植物种子萌发、营养生长、生殖器官分化、种子发育成熟过程中存在的重要调控物质。不同内源激素生理功能不同，协同调控植物生长发育。文献报道，脱落酸(ABA)通过调节细胞信号转导从而在转录水平调控下游高温胁迫诱导基因的表达，提高植物耐热性[27]。ABA与乙烯(ETH)可作为衰老信号，促进衰老提前。生长素(IAA)促进叶片生长和生理功能，后期加速衰老。玉米籽粒内源激素调节叶片光合产物向籽粒调运并积累，即灌浆速率与赤霉素(GA)、IA A、玉米素(ZR)、ABA密切相关[28]。高温胁迫改变籽粒淀粉合成酶活性。光合产物蔗糖运输到籽粒中后，在束缚态淀粉酶、可溶性淀粉合成酶、蔗糖合成酶、腺苷二磷酸葡萄糖焦磷酸化酶等一系列酶的催化下形成淀粉存储[29]。高温胁迫降低酶活性，减少淀粉含量从而引起减产。不同生长发育阶段、不同基因型玉米淀粉合成酶受高温影响不同。张保仁等[30]发现，玉米苗期高温对蔗糖磷酸合成酶、蔗糖合成酶、腺苷二磷酸葡萄糖焦磷酸化酶影响较大。水稻籽粒千粒质量下降主要是由于可溶性淀粉合成酶受到抑制[31]。高温胁迫影响产量指标特性。作物产量指标主要由穗粒数和粒质量决定。穗粒数与最大灌浆速率呈正相关，粒质量与灌浆速率最大时的生长量呈正相关。抽雄吐丝期高温胁迫降低花粉活性，雌穗小花败育增加，结实率下降[32]。不同高温胁迫持续时间对玉米的伤害程度不同，不同基因型玉米响应抽雄吐丝期高温胁迫的生理生化特性也有所不同。但目前对不同基因型玉米在不同高温胁迫时间下的耐热性分子机制尚不明确。

本研究在抽雄吐丝期对我国广泛种植的玉米品种郑单958、先玉335进行不同高温胁迫持续时间处理，利用高通量转录组测序技术及广泛靶向代谢组测序技术研究玉米响应高温胁迫的分子机理，筛选出差异表达基因、差异表达代谢物，分析其代谢通路，从而为用分子生物学技术手段提高玉米耐热性提供理论依据。

1 材料和方法

1.1 试验材料

供试材料为我国生产上广泛种植的玉米杂交种郑单958、先玉335，于2018年6月12日播种于河南省农业科学院原阳县现代农业科技示范基地。郑单958由河南省农业科学院粮食作物研究所选育，为较耐热品种。先玉335由铁岭先锋种子研究有限公司选育，并于2004年通过审定。2个品种播种均重复2次，每个品种播种20行，行长5 m，行距60 cm。试验田地力均匀，排灌方便，标准化大田管理。

1.2 高温处理

高温处理中，用长×宽×高=20 m×15 m×4 m的自制生长箱框架固定于田间，周围用透光率为95%的树脂薄膜覆盖，顶部密封90%，均匀留出10%的空隙便于气体交换。生长箱四周及中间分上部、中部悬挂温度计，以便记录处理温度。生长箱框架下部周围树脂薄膜可以拆移，当温度过高时移开薄膜降温，控制高温处理期间生长箱内气温平均高于外界田间4 ℃。生长箱内光照强度、相对湿度等生长条件基本与外界生长条件一致。

从第9片叶展开时每天8:30-17:30通过覆盖薄膜进行增温处理，至抽雄散粉7 d取第1次样。继续高温处理7 d，直至散粉结束灌浆初期取第2次样。用温湿计记录群体穗位处气温及相对湿度。正常生长条件下的郑单958和先玉335作为对照材料。取样结束后拆除高温箱使玉米材料恢复正常生长。

1.3 总RNA提取及转录组测序

高温处理7 d后玉米植株进入抽雄散粉期，取第1次样。每个样本由5株长势较为一致玉米的穗位叶混合组成，生物学重复3次。命名：HT郑单958-1、HT郑单958-2、HT郑单958-3、HT先玉335-1、HT先玉335-2、HT先玉335-3。正常生长条件下的玉米植株为对照材料，取样命名：CK郑单958-1、CK郑单958-2、CK郑单958-3、CK先玉335-1、CK先玉335-2、CK先玉335-3。高温处理14 d后，玉米植株授粉结束进入灌浆初期取第2次样，命名：HT郑单958-4、HT郑单958-5、HT郑单958-6、HT先玉335-4、HT先玉335-5、HT先玉335-6。对照材料命名：CK郑单958-4、CK郑单958-5、CK郑单958-6、CK先玉335-4、CK先玉335-5、CK先玉335-6。取样后立即放于液氮中，然后于-80 ℃超低温冰箱保存。

按照RNeasy plant mini kit(Qiagen，Germany)试剂盒操作步骤提取上述样本总RNA。用1%琼脂糖凝胶电泳检验总RNA的完整性，用Agilent bioanalyzer 2100 system(Aglilent technologies，USA)生物芯片分析检测仪检测所提RNA的质量，并用Nanodrop ND-1000 UV-Vis Spectrophotometer(Nanodrop technologie，USA)定量检测所提RNA的浓度。将检测合格的RNA 样本送至北京百迈客生物科技有限公司构建cDNA文库，然后用Illumina HiSeqTM2500测序仪进行高通量测序，结果序列为pair-end序列。

1.4 基因功能注释

通过去除接头序列、低质量序列、多N序列及长度过短序列，对测序原始序列进行质量控制得到干净序列，并计算测序碱基错误率及GC碱基含量。利用Hisat 2序列比对软件采用BWT算法将所得干净序列与参考基因组B73序列(https://www.maizegdb.org/)进行比对。利用RSeQC-2.6.3软件对所得基因进行测序饱和度、基因覆盖度及冗余序列分析后，参照Nr(NCBI non-redundant protein sequences)、Nt(NCBI non-redundant nucleotide sequence)、Pfam(Protein family)、KOG/COG(Clusters of orthologous groups of proteins)、Swiss-Prot(A manually annotated and reviewed protein sequence database)、KO(KEGG orthology database)、GO(Gene ontology)数据库，进行基因功能注释。

1.5 基因表达量计算及差异基因筛选

在转录组测序分析中用FPKM值(每100万条序列中每个基因以1 000个碱基为单位，与公布的基因组序列比对上的序列数量)衡量基因表达水平。依据基因数量并用edgeR软件计算差异表达基因，P< 0.05。显著差异表达基因筛选标准为错误发现率(FDR)< 0.05 以及基因表达量变化倍数取以2为底对数值的绝对数|log2FC|≥2。之后参照GO数据库，将差异基因按照参与的生物学过程、过程细胞组成、实现的分子功能进行基因功能分类。参照KEGG(Kyoto encyclopedia of genes and genomes)数据库分析差异表达基因参与的代谢途经。对差异表达基因的转录因子进行预测分析。

1.6 代谢物提取及检测

代谢物是生物体受遗传控制和环境影响的最终表达产物。将上述24个花粒期高温处理前后的郑单958、先玉335样本送至百迈客生物科技有限公司利用广泛靶向代谢组学技术进行代谢物定性定量分析。将叶片磨样后溶于提取液，微孔滤膜过滤后用液相色谱串联质谱(UPLC-MS/MS)检测代谢物。根据UPLC-MS/MS二级谱信息定性分析代谢物。利用三重四级杆质谱多反应检测代谢物质谱数据。

1.7 代谢组学数据分析

代谢组学数据分析主要包括差异代谢物筛选和代谢途径分析。对照公司自建代谢物数据库，用Analyst 1.6.3及Mutiaquant软件分析代谢物质谱数据。对样本主成分分析(PCA)获得样本之间总体代谢差异和组内样本之间的变异度大小，即筛选差异代谢物。采用极差法归一化处理代谢物含量，利用R软件(www.r-project.org)对代谢物在不同样本间的积累模式进行聚类分析(Hierarchical cluster analysis，HCA)。利用KEGG数据库对差异代谢物进行注释并对其代谢通路类型进行分类分析。

2 结果与分析

2.1 转录组测序统计及质量评估

转录组测序共获得214.81 Gb干净序列，24个样本的各碱基相等，没有AT/GC分离现象。识别精度99.9%以上，序列所占比例即Q30值均高于91.92%。24个样本的干净序列与玉米自交系B73网络公布基因组序列比对效率为75.34%～85.42%，详见表1。

2.2 高温胁迫相关差异表达基因检测

所得比对序列经过饱和度、覆盖度、冗余度检测后，进行可变剪接预测分析、基因结构优化分析。基于比对结果，用FPKM值衡量基因表达水平。不同处理条件下表达水平存在显著差异的基因定义为差异表达基因(DEG)。根据2组样品之间表达水平的相对高低，差异表达基因可以划分为上调基因(Up-regulated gene)和下调基因(Down-regulated gene)。通过对郑单958、先玉335高温胁迫处理7，14 d前后24个穗位叶样本转录组测序共检测到了47 026个表达发生变化的基因。郑单958抽雄期高温处理7 d后与未进行高温处理样本对比组共检测到49个表达变化水平大于2倍的差异表达基因，其中24个为上调表达基因，25个为下调表达基因。继续高温处理14 d后与未进行高温处理样本对比组共检测到306个差异表达基因，其中130个为上调表达基因，176个为下调表达基因。高温处理7 d样本与高温处理14 d样本对比组共检测到1 462个差异表达基因，其中647个为上调表达基因，815个为下调表达基因。先玉335高温处理7 d后与未进行高温处理样本对比组共检测到381个差异表达基因，其中164个为上调表达基因，217个为下调表达基因。继续高温处理14 d后与未进行高温处理样本对比组检测到299个差异表达基因，其中226个为上调表达基因，73个为下调表达基因。高温处理7 d样本与高温处理14 d样本对比组检测到2 481个差异表达基因，其中1 275个为上调表达基因，1 206个为下调表达基因。郑单958和先玉335品种间差异表达基因比较多，高温处理7 d后检测到6 646个差异表达基因，其中3 253个为上调表达基因，3 393个为下调表达基因。高温处理14 d后检测到5 958个差异表达基因，其中3 110个为上调表达基因，2 848个为下调表达基因(表2)。

表1 24个玉米花粒期高温胁迫样本转录组测序数据

表2 花粒期高温胁迫处理不同对比分组中差异表达基因数目统计

2.3 高温胁迫相关差异表达基因分析

参照KEGG pathway系统基因组信息数据库，对检测到的差异表达基因代谢通路进行富集分析。KEGG代谢通路数据库包括碳水化合物代谢、核苷酸及氨基酸代谢、有机物生物合成及降解等多种代谢途经。郑单958 高温胁迫7 d后检测到的49个差异表达基因主要分布于淀粉和蔗糖代谢、玉米素生物合成、牛磺酸及亚磺酸代谢、丁酸甲酯代谢、β-丙氨酸代谢等代谢途经中。高温胁迫14 d后差异表达基因主要分布于苯丙素生物合成、甘油磷脂代谢、淀粉及蔗糖代谢、硫胺素代谢、核黄素代谢、氮代谢等代谢途经中。从图1可以看出，不同高温胁迫对比组中检测到了更多的差异表达基因，主要集于内质网蛋白质加工、谷胱甘肽代谢、苯并噁嗪类化合物生物合成、其他类型氧化多糖生物合成、植物昼夜节奏、氮代谢、酪氨酸代谢等代谢途经。先玉335高温胁迫7 d后差异表达基因主要分布于内质网蛋白质加工、淀粉及蔗糖代谢、戊糖与葡萄糖醛酸酯相互转化、植物-病原体相互作用、玉米素生物合成、苯并噁嗪类化合物生物合成、苯丙素生物合成、植物激素信号传导等代谢途经中。高温胁迫14 d后差异表达基因主要分布于苯丙素生物合成、植物昼夜节奏、植物激素信号传导、淀粉及蔗糖代谢、光合天线蛋白、硫胺素代谢、亚油酸代谢、酪蛋白氨基酸代谢、甘氨酸，丝氨酸及苏氨酸代谢等代谢途径中。先玉335不同高温胁迫时间样品对比组检测到了更多的差异表达基因，如图2所示，主要富集于内质网蛋白质加工代谢、植物激素信号传导、植物-病原体相互作用、类胡萝卜素生物合成、谷胱甘肽代谢、植物昼夜节奏、精氨酸及脯氨酸代谢、其他类型氧化多糖生物合成等代谢途经。郑单958与先玉335品种间高温胁迫7 d对比组中检测到的差异表达基因主要集中于谷胱甘肽代谢、氮代谢、脂肪酸生物合成、鞘脂类化合物代谢、黄酮类化合物生物合成、黄酮及黄酮醇生物合成、类单萜生物合成、维生素B6代谢等代谢途经。郑单958与先玉335品种间高温胁迫14 d所检测到的差异表达基因主要集中于谷胱甘肽代谢、丙酮酸代谢、半胱氨酸及蛋氨酸代谢、鞘脂类化合物代谢、氮代谢、精氨酸及脯氨酸代谢、硫代谢、苯基丙氨酸代谢、醚酯代谢、花生四烯酸代谢等代谢途经中。

每一个圆点代表一个KEGG通路；纵坐标代表通路名称，横坐标为富集因子，即差异基因中注释到某通路的基因比例与所有基因中注释到该通路的基因比例的比值；富集因子越大，表示差异表达基因在该通路中的富集水平越显著；圆点颜色代表q值，从上到下代表差异表达基因在该通路中的富集显著性越可靠；圆点大小代表通路中富集的差异表达基因数目，圆点越大基因越多。图2同。

Each dot represents a KEGG path； The ordinate represents the name of the pathway，and the abscissa is the enrichment factor;The greater the enrichment factor，the more significant the enrichment level of differentially expressed genes in this pathway;The dot color represents the value ofqvalue. The lower the dot，the more reliable the enrichment significance of differentially expressed genes in this pathway;the size of dots represents the number of differentially expressed genes enriched in the pathway.The larger the dots，the more genes there are.The same as Fig.2.

图1 HT郑单958-1-2-3 vs HT郑单958-4-5-6对比组差异表达基因KEGG通路富集散点图

Fig.1 Differentially expressed gene KEGG pathway enrichment scatter map ofHT Zhengdan 958-1-2-3 vs HT Zhengdan 958-4-5-6 group

2.4 高温胁迫相关转录因子分析

通过对转录组测序后转录因子分析，共获得67个不同转录因子类型。其中主要有AP2/ERF、AGC、MYB、Aux/IAA、B3、C2C2、CAMK、CMGC、AGARP、HSF、WRKY、SWI、TKL、BHLH、C3H等15类参与植物生长发育调控、抗逆胁迫相关重要转录因子。Aux/IAA参与生长素调节，在抗旱、耐盐中起重要作用。AP2/ERF转录因子广泛存在植物体内，参与生长发育及非生物胁迫相关基因的表达调控。热激蛋白转录因子HSF调控下游基因直接参与植物高温胁迫应答网络。

图2 HT先玉335-1-2-3 vs HT先玉335-4-5-6对比组差异表达基因KEGG通路富集散点图

2.5 高温胁迫相关差异代谢物筛选

高温胁迫后玉米植株内代谢产物产生不同的应答机制。代谢组学分析可以从高温胁迫样本中筛选出具有生物学意义和统计显著差异的代谢物，进而阐明应答高温胁迫的代谢过程和变化机理。本试验共检测到654个代谢物，采用正交偏最小二乘法-判别分析统计方法分析不同代谢物在试验样本组间差异。最终以P< 0.05且|log2FC|≥2为标准筛选差异代谢物(DEM)，如表3所示。郑单958高温胁迫7 d之后与正常生长组比较共检测出28个差异代谢物，注释到5个代谢通路中，其中9个差异代谢物呈上调表达，19个下调表达。郑单958高温胁迫14 d之后与正常生长组比较共检测出54个差异代谢物，注释到13个代谢通路中，其中9个差异代谢物呈上调表达，45个下调表达。先玉335高温胁迫7 d之后与正常生长组比较共检测出98个差异代谢物，注释到24个代谢通路中，其中43个差异代谢物呈上调表达，55个下调表达。先玉335高温胁迫14 d之后与正常生长植株对比组共检测出38个差异代谢物，注释到13个代谢通路中，其中14个差异代谢物呈上调表达，24个下调表达。郑单958高温胁迫7 d与先玉335高温胁迫7 d对比组共检测出144个差异代谢物，注释到36个代谢通路中，其中81个差异代谢物呈上调表达，63个下调表达。郑单958高温胁迫14 d与先玉335高温胁迫14 d对比组共检测出158个差异代谢物，注释到40个代谢通路中，其中81个差异代谢物呈上调表达，77个下调表达。

2.6 不同玉米品种高温胁迫相关差异表达基因分析

对比分析郑单958、先玉335品种间不同高温胁迫时间对比组，检测到差异表达基因。CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组与CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组差异表达基因韦恩图如图3所示，检测到8个共同差异表达基因，均在郑单958和先玉335高温胁迫7 d后的叶片中表达。分别为：Zm00001d015825，编码玉米果胶裂解酶超级家族蛋白；Zm00001d018279，编码玉米LOC100285574非特性蛋白前体；Zm00001d023387，编码玉米LOC100278780非特性蛋白；Zm00001d032046，编码玉米细胞分裂素脱氢酶10；Zm00001d040171，编码玉米LOC100192465蛋白；Zm00001d044728，编码玉米热激蛋白；Zm00001d047146，编码玉米膨胀蛋白-B1前体；Zm00001d048592，编码玉米核酮糖二磷酸羧化酶激酶1。CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组中检测到41个特异差异表达基因。CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组中检测到373个特异差异表达基因，远远多于郑单958高温胁迫7 d对比组。推测高温胁迫7 d后更多的差异表达基因参与高温胁迫响应网络。参照KEGG pathway数据库对比分析发现，编码玉米果胶裂解酶超级家族蛋白的Zm00001d015825及编码玉米LOC100285574非特性蛋白前体的Zm00001d0182792个共同差异表达基因富集在淀粉和蔗糖代谢途径中，在郑单958和先玉335对比组中均呈下调表达。

表3 花粒期高温胁迫处理不同对比分组中差异代谢物数目

G0.CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组；G4.CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组；G1.CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组；G5.CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组；G6.HT郑单958-1-2-3 vs HT郑单958-4-5-6对比组；G8.HT先玉335-1-2-3 vs HT先玉335-4-5-6对比组。

G0.CK Zhengdan 958-1-2-3 vs HT Zhengdan 958-1-2-3 contrast group；G4.CK Xianyu 335-1-2-3 vs HT Xianyu 335-1-2-3 contrast group；G1.CK Zhengdan 958-4-5-6 vs HT Zhengdan 958-4-5-6 contrast group；G5.CK Xianyu 335-4-5-6 vs HT Xianyu 335-4-5-6 contrast group；G6.HT Zhengdan 958-1-2-3 vs HT Zhengdan 958-4-5-6 contrast group；G8.HT Xianyu 335-1-2-3 vs HT Xianyu 335-4-5-6 contrast group.

图3 不同品种对比组中差异表达基因韦恩图

Fig.3 Venn diagram of differentially expressed genes in different cultivar groups

CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组与CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组中，检测到14个共同差异表达基因，均在郑单958和先玉335高温胁迫14 d后的叶片中表达。高温胁迫14 d比高温胁迫7 d对比组中共同差异基因略多，推测是由于更多的基因参与响应花粒期玉米高温胁迫引起的。分别为Zm00001d010410，编码玉米水通道TIP4-2 相似蛋白；Zm00001d010662，编码玉米SNF1 相关蛋白激酶调控亚蛋白β-1亚型；Zm00001d011183，编码玉米叶绿体前体中硫胺素噻唑合成酶1；Zm00001d013781，编码玉米叶绿体光吸收叶绿素荧光增强蛋白；Zm00001d014680，编码玉米阳离子转运体HKT9；Zm00001d015434，编码玉米B-环锌指家族蛋白；Zm00001d018234，编码玉米非特性LOC100275888蛋白；Zm00001d028154，编码玉米转移酶亚型Ⅱ；Zm00001d033869，编码玉米三角状五肽重复结构域蛋白At3g49730；Zm00001d038746，编码玉米热激转录受体A-4d 亚体 X2；Zm00001d041776，编码玉米B-葡糖苷酶；Zm00001d043050，编码玉米非特性LOC111347718 蛋白；Zm00001d049006，编码玉米PHD锌指蛋白；Zm00001d053042，编码玉米非特性LOC100277850 蛋白。CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组、CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组分别有292，285个特异的差异表达基因。郑单958高温胁迫14 d样品比先玉335高温胁迫14 d样本差异表达基因略多，推测有更多的差异表达基因参与郑单958高温胁迫响应机制。通过KEGG富集分析，共同差异表达基因主要分布于硫胺素代谢、硫代谢、氰基氨基酸代谢、半胱氨酸及蛋氨酸代谢、苯丙素生物合成、淀粉及蔗糖代谢、碳代谢8个代谢途经中。Zm00001d011183富集在硫胺素代谢途径中，在郑单958和先玉335对比组中均呈上调表达。

HT郑单958-1-2-3 vs HT郑单958-4-5-6对比组与HT先玉335-1-2-3 vs HT先玉335-4-5-6对比组中检测到756个共同差异表达基因。HT郑单958-1-2-3 vs HT郑单958-4-5-6对比组、HT先玉335-1-2-3 vs HT先玉335-4-5-6对比组中分别检测到706，1 725个特异差异表达基因，远远多于高温胁迫7 d、高温胁迫14 d检测到的共同差异表达基因和特异差异表达基因。推测随着高温胁迫时间延长，玉米植株体内更多的差异表达基因参与响应高温胁迫。通过KEGG pathway数据库分析，共同差异表达基因主要富集于内质网蛋白质加工过程、谷胱甘肽代谢、其他类型O-多糖生物合成、氨基多糖降解、苯并噁嗪类生物合成、其他聚糖降解、植物昼夜节律等20个代谢途经中。

2.7 高温胁迫时间对玉米响应高温相关差异表达基因的影响

从图4可以看出，CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组与CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组中没有检测到共同差异表达基因。CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组、CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组分别检测到49，306个特异差异表达基因。郑单958随着高温胁迫的延长，更多的差异表达基因参与了响应高温胁迫。郑单958不同高温胁迫时间对比组中没检测到共同差异表达基因，可能是由于高温胁迫7 d对比组、高温胁迫14 d对比组中差异表达基因均比较少造成的。

CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组与CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组中检测到12个共同差异表达基因。这12个共同差异表达基因分别为：Zm00001d005750，编码玉米转录因子bHLH69；Zm00001d011486，编码玉米CTP合成酶异构体X3；Zm00001d012378，编码玉米非特性LOC100274367蛋白；Zm00001d017422，编码玉米同源异型-亮氨酸框拉链蛋白ATHB-6；Zm00001d017666，编码玉米NRT1/PTR家族蛋白异构体X1；Zm00001d019399，编码玉米类PVR3蛋白前体；Zm00001d019565，编码玉米水通道蛋白PIP2-6；Zm00001d022069，编码玉米忍耐蛋白5；Zm00001d025037，编码玉米UDP-糖基转移酶79；Zm00001d033543，编码玉米跨膜HPP 家族蛋白；Zm00001d049789，编码玉米激素应激家族蛋白；Zm00001d051362，编码玉米水通道TIP2-1蛋白。12个共同差异表达基因主要富集于吡啶代谢通路。KEGG pathway数据库对比分析注释共同差异表达基因Zm00001d011486富集在吡啶代谢途径中，2个对比组中均呈下调表达。CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组、CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组分别有369，287个特异差异表达基因。先玉335高温胁迫14 d后检测到的差异表达基因比高温胁迫7 d少，推测先玉335抗高温调控网络主要发生在高温胁迫初期，与生产上先玉335比郑单958抗高温性差一致。

HT郑单958-1-2-3 vs HT先玉335-1-2-3对比组与HT郑单958-4-5-6 vs HT先玉335-4-5-6对比组中检测到3 767个共同差异表达基因。在HT郑单958-1-2-3 vs HT先玉335-1-2-3对比组检测到了2 879个特异差异表达基因，在HT郑单958-4-5-6 vs HT先玉335-4-5-6对比组检测到2 191个特异差异表达基因。推测这是由于高温胁迫时间差异及玉米品种差异共同造成的。通过KEGG pathway数据库富集分析，差异表达基因主要分布于谷胱甘肽代谢，鞘脂类代谢，黄酮、黄酮醇类化合物生物合成，氮代谢，丙酮酸代谢，硫代谢，花生四烯酸代谢，光合作用天线蛋白，类黄酮醇生物合成，缬氨酸、亮氨酸及异亮氨酸生物合成等20个代谢途径中。KEGG pathway数据库对比分析发现，有24个共同差异基因富集在谷胱甘肽代谢途径中，17个共同差异基因富集在丙酮酸代谢途径中。

G0.CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组；G1.CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组；G4.CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组；G5.CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组；G2.HT郑单958-1-2-3 vs HT先玉335-1-2-3对比组；G3.HT郑单958-4-5-6 vs HT先玉335-4-5-6对比组。

G0.CK Zhengdan 958-1-2-3 vs HT Zhengdan 958-1-2-3 contrast group；G1.CK Zhengdan 958-4-5-6 vs HT Zhengdan 958-4-5-6 contrast group；G4.CK Xianyu 335-1-2-3 vs HT Xianyu 335-1-2-3 contrast group；G5.CK Xianyu 335-4-5-6 vs HT Xianyu 335-4-5-6 contrast group；G2.HT Zhengdan 958-1-2-3 vs HT Xianyu 335-1-2-3 contrast group；G3.HT Zhengdan 958-4-5-6 vs HT Xianyu 335-4-5-6 contrast group.

图4 不同时间对比组中差异表达基因韦恩图

Fig.4 Venn diagram of differentially expressed genes in different time groups

2.8 不同玉米品种间高温胁迫相关差异代谢物表达分析

CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组与CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组共筛选到13个共同差异代谢物，均在2个对比组中表达，分别为：pmb0645，6-C-己糖基-橙皮素邻己糖苷；pmb3074，3-O-P-香豆蔻酰奎宁酸；pmb3093，6，7-二羟基香豆素7-O-奎宁酸； pmb3094，阿舒利汀邻喹酰阿舒利汀邻喹酸；pmb3121，LysoPE 18∶1；pmb3132，LysoPE 18∶1异构体；pme0199，槲皮素；pme0399，绿原酸(3-O-咖啡酰奎宁酸)；pme1539，异鼠李素3-O-新橙皮苷；pme1541，阿卡西汀；pme1816，新绿原酸(5-O-咖啡酰奎宁酸)；pme2901，1-O-咖啡酰奎宁酸；pme3514，桑色素。CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组、CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组中分别检测到15，85个特异差异表达代谢物。郑单958高温胁迫处理7 d后叶片中的差异表达代谢物明显少于先玉335。通过与KEGG pathway数据库对比分析，检测的所有差异表达代谢物富集到20个主要的代谢通路上，分别为：类黄酮生物合成，氧化磷酸化，光合作用，黄酮、黄酮醇生物合成，烟酸、烟酰胺代谢，苯丙氨酸代谢，苯丙素类生物合成，丙氨酸、天门冬氨酸及谷氨酸代谢，丁酸盐代谢，酪氨酸代谢，柠檬酸循环，碳青霉烯类生物合成，泛醌及其他类萜醌类生物合成，维生素B6代谢，芪类、二芳基庚酸类及姜辣素生物合成，不饱和脂肪酸生物合成，叶酸生物合成，α-亚麻酸代谢，核黄素代谢，丙酮酸代谢。

CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组与CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组检测到12个共同差异代谢物，分别为：pmb0130，N-乙酰色胺；pmb0490，N-P-腐胺；pmb0505，N′-P-腐胺；pmb0774，N-羟色胺；pmb0883，LysoPE 18∶0；pmb2561，N-乙酰蛋氨酸；pme1841，1，5-二氨基戊烷；pme2024，5-羟基色胺；pme2651，NADP；pme3005，N，N-二甲基甲酰胺；pme3038，5-奥托洛林；pme3104，1-甲基腺嘌呤。CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组、CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组中分别有42，26个特异差异表达代谢物。郑单958高温胁迫处理14 d后叶片中的差异表达代谢物多于先玉335，这与郑单958高温胁迫14 d差异表达基因多于先玉335的检测结果一致。通过KEGG pathway数据库注释分析，差异代谢物富集在代谢途径，次生代谢产物生物合成，2-氧羧酸代谢，嘌呤代谢，异黄酮生物合成，光合作物，谷胱甘肽代谢，赖氨酸降解，芪类、二芳基庚酸类及姜辣素生物合成，丙酮酸代谢，玉米素生物合成，类黄酮生物合成，色氨酸代谢，烟酸及烟酰胺代谢，黄酮、黄酮醇生物合成等20个代谢通路上。

HT郑单958-1-2-3 vs HT郑单958-4-5-6对比组与HT先玉335-1-2-3 vs HT先玉335-4-5-6对比组共筛选到16个共同差异表达代谢物，均在2个对比组中表达(表4)。HT郑单958-1-2-3 vs HT郑单958-4-5-6对比组、HT先玉335-1-2-3 vs HT先玉335-4-5-6对比组分别检测到33，60个特异差异表达代谢物。先玉335不同高温胁迫时间后检测到的差异表达代谢物多于郑单958不同高温胁迫时间样本。通过KEGG pathway数据库注释，差异表达代谢物富集到类黄酮代谢，氰氨基酸生物合成，半胱氨酸及蛋氨酸代谢，精氨酸生物合成，异黄酮生物合成，单环丙酰胺生物合成，玉米素生物合成，不饱和脂肪酸生物合成，丙氨酸、天门冬氨酸及谷氨酸代谢等20个代谢通路上。

2.9 高温胁迫时间对差异表达代谢物的影响

从图5可以看出，CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组与CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组共筛选到11个共同差异表达代谢物，均在2个对比组中表达。分别为：pmb0896，N′-阿魏酰腐胺；pmb0907，N-阿魏酰腐胺；pmb3074，3-O-P-香豆蔻酰奎宁酸；pmb3093，6，7-二羟基香豆素7-O-奎宁酸；pmb3094，阿舒利汀邻喹酰阿舒利汀邻喹酸；pme0199，槲皮素；pme0399，绿原酸(3-O-咖啡酰奎宁酸)；pme1683，D-泛酸；pme1816，新绿原酸(5-O-咖啡酰奎宁酸)；pme2901，1-O-咖啡酰奎宁酸；pme3514，桑色素。CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组、CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组中分别检测到17，43个特异差异表达代谢物。随着高温胁迫时间的延长，郑单958植株内差异表达代谢物也随着增加。通过KEGG pathway数据库对比分析，差异代谢物富集到代谢途径，黄酮及黄酮醇生物合成，甘氨酸、丝氨酸及苏氨酸代谢，次生代谢产物生物合成，氨基酸生物合成，氨酰-tRNA生物合成，类黄酮生物合成，精氨酸及脯氨酸代谢，碳代谢，半胱氨酸及蛋氨酸代谢，牛磺酸及亚磺酸代谢11个代谢通路。

CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组与CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组共筛选到11个共同差异代谢物，均在2个对比组中表达，分别为：pma0253，O-甲基金圣草黄素 5-O-己糖苷；pma6638，O-甲基金圣草黄素7-O-己糖苷；pmb2999，金圣草黄素 5-O-己糖苷；pmb3114，表儿茶素-表阿夫儿茶精；pme0088，木犀草素；pme0394，2-甲氧基苯甲酸；pme1622，山柰酚3-O-葡萄糖苷；pme2651，NADP；pme3104，1-甲基腺嘌呤；pme3276，2′-羟基染料木素；pme3367，6-甲基香豆素。CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组、CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组中分别检测到87，27个特异差异表达代谢物。先玉335中随着高温胁迫处理时间的延长，差异表达代谢物反而减少，这与郑单958中差异表达代谢物检测结果相反。生产上先玉335比郑单958耐高温性差[33]，推测正是由于先玉335中随着高温胁迫时间延长响应高温胁迫的代谢物反而减少的缘故。通过KEGG pathway数据库注释，差异代谢物富集到类黄酮生物合成，氧化磷酸化，黄酮及黄酮醇生物合成，烟酸及烟酰胺代谢，碳青霉烯代谢，柠檬酸循环，苯丙素生物合成，丙氨酸、天门冬氨酸及谷氨酸代谢等20个代谢通路上。

表4 HT郑单958-1-2-3 vs HT郑单958-4-5-6对比组与HT先玉335-1-2-3 vs HT先玉335-4-5-6对比组中共同差异表达代谢物

G4.CK郑单958-1-2-3 vs HT郑单958-1-2-3对比组；G5.CK郑单958-4-5-6 vs HT郑单958-4-5-6对比组；G2.CK先玉335-1-2-3 vs HT先玉335-1-2-3对比组；G3.CK先玉335-4-5-6 vs HT先玉335-4-5-6对比组；G9.HT先玉335-1-2-3 vs HT郑单958-1-2-3对比组；G11.HT先玉335-4-5-6 vs HT郑单958-4-5-6对比组。

G4.CK Zhengdan 958-1-2-3 vs HT Zhengdan 958-1-2-3 contrast group；G5.CK Zhengdan 958-4-5-6 vs HT Zhengdan 958-4-5-6 contrast group；G2.CK Xianyu 335-1-2-3 vs HT Xianyu 335-1-2-3 contrast group；G3.CK Xianyu 335-4-5-6 vs HT Xianyu 335-4-5-6 contrast group；G9.HT Xianyu 335-1-2-3 vs HT Zhengdan 958-1-2-3 contrast group；G11.HT Xianyu 335-4-5-6 vs HT Zhengdan 958-4-5-6 contrast group.

图5 不同对比组间差异表达代谢物韦恩图

Fig.5 Venn diagram of differentially expressed metabolites in different groups

HT先玉335-1-2-3 vs HT郑单958-1-2-3对比组与HT先玉335-4-5-6 vs HT郑单958-4-5-6对比组共筛选到95个共同差异代谢物，均在2个对比组中表达。HT先玉335-1-2-3 vs HT郑单958-1-2-3对比组、HT先玉335-4-5-6 vs HT郑单958-4-5-6对比组中分别检测到49，63个特异差异表达代谢物。郑单958、先玉335品种间差异及不同高温胁迫时间共同引起了更多代谢物表达发生变化[34]，植物响应高温胁迫机制与基因型、胁迫时间密切相关。通过KEGG pathway数据库对比分析，差异代谢物富集到嘌呤代谢，精氨酸生物合成，黄酮及黄酮醇生物合成，色氨酸代谢，牛磺酸及亚磺酸代谢，烟酸及烟酰胺代谢，氨酰tRNA生物合成，甘氨酸、丝氨酸及苏氨酸代谢，泛酸盐及辅酶A生物合成等20个代谢通路上。

3 结论与讨论

通过对郑单958、先玉335高温胁迫处理7，14 d前后24个穗位叶样本转录组测序，共检测到了47 026个表达发生变化的基因。进一步筛选发现，郑单958高温胁迫处理7 d后有49个差异表达基因，24个呈上调表达。高温胁迫14 d后检测到306个差异表达基因，130个呈上调表达。随着高温胁迫时间的延长，郑单958中差异表达基因增多。先玉335高温胁迫处理7 d后有381个差异表达基因，164个呈上调表达。高温胁迫14 d后检测到299个差异表达基因，226个呈上调表达。可见不同基因型玉米响应高温胁迫的内在机制不同。差异表达基因韦恩图表明，郑单958高温胁迫7 d与先玉335高温胁迫7 d对比组中共检测到8个共同差异表达基因。差异表达基因主要富集于淀粉及蔗糖代谢、硫胺素代谢等通路中。郑单958与先玉335高温胁迫14 d对比组中检测到14个共同差异表达基因。主要富集于内质网内蛋白质加工过程、淀粉及蔗糖代谢、戊糖及葡萄糖醛酸酯相互转换、苯丙素生物合成、植物昼夜节律、光合天线蛋白、硫胺素代谢等通路中。郑单958 vs先玉335高温胁迫7 d 与郑单958 vs先玉335高温胁迫14 d对比组中共检测到3 767个共同差异表达基因。主要富集于谷胱甘肽代谢、鞘脂类化合物生物合成、氮代谢、类黄酮生物合成、硫代谢、丙酮酸代谢、苯丙氨酸代谢等代谢通路中。可见基因型及高温胁迫时间对差异表达基因有着较大的影响。

通过对郑单958、先玉335高温胁迫处理7，14 d前后24个穗位叶样本代谢组测序共检测到了654个代谢物。进一步筛选发现郑单958高温胁迫处理7 d后有28个差异表达代谢物，9个呈上调表达。高温胁迫14 d后检测到54个差异表达代谢物，9个呈上调表达。随着高温胁迫时间的延长，郑单958中差异表达代谢物增多。先玉335高温胁迫处理7 d后有98个差异表达代谢物，43个呈上调表达。高温胁迫14 d后检测到38个差异表达代谢物，14个呈上调表达。随着高温胁迫时间的延长，先玉335中差异表达代谢物反而减少。差异表达代谢物韦恩图表明，郑单958高温胁迫7 d前后与高温胁迫14 d前后对比组中共检测到11个共同差异表达代谢物。主要富集于类黄酮代谢、黄酮及黄酮醇生物合成、烟酸及烟酰胺代谢、苯丙素生物合成、氧化磷酸化等代谢通路中。先玉335高温胁迫7 d前后与先玉335高温胁迫14 d前后对比组中共检测到11个共同差异表达代谢物。主要富集于代谢路径、次生代谢产物生物合成、嘌呤代谢、异黄酮生物合成、谷胱甘肽代谢等代谢通路中。郑单958 vs先玉335高温胁迫7 d 与郑单958 vs先玉335高温胁迫14 d对比组中共检测到95个共同差异表达代谢物。主要富集于嘌呤代谢、精氨酸生物合成、氨酰tRNA-生物合成、色氨酸代谢、黄酮及黄酮醇生物合成、烟酸及烟酰胺生物合成等代谢通路中。可见基因型及高温胁迫时间共同影响玉米响应高温的调控网络。