轉錄組及代謝組聯合解析玉米響應花粒期高溫脅迫機制

李 川，喬江方，黃 璐，張美微，張盼盼，牛 軍，劉京寶

(河南省農業科學院 糧食作物研究所，河南 鄭州 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-生物合成、色氨酸代謝、黃酮及黃酮醇生物合成、煙酸及煙酰胺生物合成等代謝通路中。可見基因型及高溫脅迫時間共同影響玉米響應高溫的調控網絡。