高溫脅迫下不同玉米材料間的基因差異表達分析

曹丹，張紅梅，楊海鵬，趙歡，劉小紅

(1.西華師范大學生命科學學院，四川 南充 637000；2.山西省農業科學院玉米研究所，山西 忻州 034000)

玉米(Zeamays)是世界上最重要的3大谷類作物之一，與水稻和小麥并列，在世界農業中占有重要地位[1].高溫是限制全球玉米產量的最關鍵的非生物脅迫因素之一，世界上許多國家的玉米生產都遭受高溫危害[2-5].在中國，黃淮海地區是我國夏玉米集中產區，幾乎每年都會遭遇高溫危害.近年來，伴隨著全球氣候變暖，干旱和高溫協同脅迫已經成為該地區玉米生產中頻繁遇到的一種自然災害，且這一危害幾乎發生在玉米整個生育期，嚴重影響著玉米產量和品質[6].

高溫脅迫可引起植物的生理、分子和生化變化，干擾細胞和整個植物的生長發育過程，從而對作物的生產產生不利影響[7].對玉米來說，除了在其他生長階段可能遭遇高溫脅迫外，在生殖生長及受精結實階段特別容易遇到高溫天氣，導致玉米減產[8].高溫對花粉(雄配子)發育的影響要大于花絲(雌配子)[9].在雌雄穗分化至籽粒成熟階段如果遭遇持續的高溫脅迫，可能造成敗育小花數量增加、花粉結構破壞、雌穗結實率降低，最終導致減產.研究發現，不同玉米材料耐高溫脅迫能力存在差異[10-11]，目前在國內耐高溫的玉米材料較少，僅有鄭單958和中地88等較少的耐高溫玉米品種被推廣應用.因此，以不同耐高溫脅迫能力的玉米花粉作為材料，研究在高溫環境下基因的差異表達具有重要意義.

基于此，本研究以耐高溫脅迫的中地88和高溫敏感的先玉335兩個玉米品種作為試驗材料，在高溫脅迫環境下取其花粉作轉錄組測序分析，了解基因的差異表達情況，并進一步對差異表達基因作功能注釋和代謝通路分析，以期為進一步闡明玉米耐高溫脅迫的分子遺傳機制奠定理論基礎，為克隆應用玉米耐高溫脅迫基因的分子育種提供參考.

1 材料與方法

1.1 試驗材料

本試驗選用了兩個玉米材料，一個是耐高溫脅迫的中地88，另一個是對高溫敏感的先玉335，兩個玉米材料均來自山西省農業科學院玉米研究所.

1.2 試驗方法

2019年將中地88和先玉335兩個玉米品種種植于山西省農業科學院玉米研究所實驗大棚，每個材料種植3行，行長3 m，每行10株.從播種到開始散粉前第5天，未作增溫處理，玉米生長環境包括溫度、濕度和光照等條件與田間自然生長的玉米材料一樣.在開始散粉前的第4天作增溫處理，但濕度和光照等條件與室外保持一致，在盛花期頭一天下午對雄穗進行套袋，第2天12時30分時監測到雄穗所處位置的溫度達到42 ℃，通過大棚通風手段維持該溫度1 h，在13時30分，取兩個材料的花粉，各取4個重復，每個重復取約400 mg花粉裝進2 mL離心管，立即液氮速凍，該過程在實驗大棚內完成.然后在實驗室提取花粉總RNA，對符合要求的總RNA每個材料各選3份(作為3個生物學重復)，中地88作為處理，用T1、T2和T3表示；先玉335作為對照，用CK1、CK2和CK3表示；再進一步完成mRNA純化、反轉錄、建庫及測序等工作(該工作由上海派森諾生物科技股份有限公司完成).

1.3 數據統計分析

1.3.1 參考基因組選擇 試驗結果分析選取玉米B73作為參考基因組，該參考物種在數據庫中注釋的基因數目及注釋比例為：NCBI_GI：16 607(41.94%)；KEGG：6 715(16.96%)；GO：30 113(76.06%)；EC：2 921(7.37%)；UniProtID：16 263(41.07%)；Ensembl：39 591(100%).

1.3.2 序列比對分析 對試驗返回的結果，使用TopHaT2的升級版HISAT2(http://ccb.jhu.edu/software/hisaT2/index.shtml) 軟件分析，將過濾后的Reads比對到參考基因組上，并對獲得的基因在基因組上的區域分布進行統計.

1.3.3 基因表達量分析 使用HTSeq統計比對到每一個基因上的Read Count值作為基因的原始表達量.采用FPKM(Fragments Per Kilobase Million)對表達量進行標準化，再繪制表達基因的FPKM密度分布的小提琴圖.

1.3.4 基因差異表達分析 根據基因在兩種玉米材料中的表達量，以先玉335材料為對照，以中地88材料為處理，分析上調基因、下調基因及無顯著差異表達基因的數量，并繪制差異表達基因的火山圖和在不同染色上的基因組圈圖.

對差異表達基因的生物學功能進行GO(gene ontology)富集及功能注釋，對差異表達基因的GO富集分析結果，按照分子功能(molecular function，MF)、生物過程(biological process，BP) 和細胞組分(cellular component，CC) 進行GO 分類，每個GO 分類中挑選P-value最小，即富集最顯著的前10個GO term條目進行展示.此外，還對差異表達基因作KEGG(kyoto encyclopedia of genes and genomes)富集.根據差異表達基因的KEGG富集分析結果，挑選P-value最小即富集最顯著的前20個Pathway進行展示.

2 結果與分析

2.1 Reads與參考基因組的比對

2.1.1 基本信息比對 使用TopHaT2的升級版HISAT2軟件將過濾后的 Reads 比對到參考基因組上.序列比對的基本信息統計結果見表1.對T1、T2、T3、CK1、CK2和CK3面言，比對上參考基因組的序列總數占比對的序列總數的百分比分別為92.50%、92.83%、92.24%、90.62%、91.25%和92.54%，平均值為91.96%；比對到多個位置的序列總數占比對上參考基因組的序列總數的百分比分別為5.88%、5.87%、6.29%、6.69%、6.47%和6.69%，平均值為6.34%；只比對到一個位置的序列總數占比對上參考基因組的序列總數分別為94.12%、94.13%、93.71%、93.31%、93.53%和93.31%，平均值為93.66%.

表1 Reads與參考基因組間的基本信息比對Table 1 Basic information comparison between reads and reference genome

2.1.2 比對區域分布 將比對到基因組上的Reads分布情況進行統計，結果見表2.對T1、T2、T3、CK1、CK2和CK3面言，比對到基因區域的Reads總數占比對事件發生的總次數的百分比分別為94.24%、94.07%、93.53%、93.76%、94.10%和94.07%，平均值為93.96%；比對到基因間區的Reads總數占比對事件發生的總次數的百分比分別為5.76%、5.93%、6.47%、6.24%、5.90%和5.93%，平均值為6.04%；比對到外顯子區域的Reads總數占比對到基因區域的 Reads 總數的百分比分別為99.57%、99.59%、99.56%、99.45%、99.50%和99.45%，平均值為99.52%.

表2 比對到參考基因組上的Reads分布Table 2 Distribution of the Reads mapped to reference genome

2.2 基因表達的結果與分析

2.2.1 基因表達量 采用FPKM對表達量進行標準化后，將表達量分為不同的區間，對各樣本在不同表達量區間內基因的數目進行統計，結果見表3.表達量值在0～0.01區間的基因數量最多，所有6個樣本都超過了20 000個，平均為21 228個，占總量的60.96%；隨表達量范圍值的提高，對應基因數量逐漸減少，在>1000表達量區間所有6個樣本的基因數量都低于160個，平均為155.5個，僅占總量的0.45%.由此表得出，低表達的基因數量較多，而高表達的基因數量較少.

表3 不同表達量區間的基因數量統計Table 3 Statistics of the number for the genes in different expression quantity intervals

2.2.2 表達基因的FPKM密度分布 基因在各個樣品中的表達量利用小提琴圖進行統計分析，結果如圖1所示.從圖1可以看出，這6個樣本具有相似的基因表達模式，即中等表達的基因占絕大多數，而低表達和高表達的基因僅占一小部分.

2.3 基因差異表達結果與分析

2.3.1 差異表達基因數量 中地88(Case:Traitment)相對于先玉335(Control：CK)材料而言，差異表達基因共計有1 822個，相比于Control，上調表達基因有857個，下調表達基因有965個.另還有22 731個基因在這兩個玉米材料中表達差異不顯著.繪制的差異表達基因的火山圖見圖2.從圖2可看出，左側(藍色)為Case相比于Control下調基因，右側(紅色)為Case相比于Control上調基因，左右差異基因分布大致對稱.

2.3.2 差異表達基因的染色體分布 對在這兩個玉米材料中表達的基因在染色體上的分布作基因組圈圖，結果如圖3所示.從上調基因、下調基因及無顯著差異表達基因在整個玉米10條染色體上的分布來看，分布較均勻.

x為FPKM值；圖中盒型中間的橫線是中位數，盒型的上下邊緣為75%，上下限為90%，外部形狀為核密度估計.x means FPKM value；The horizontal line in the middle of the box is the median，the upper and lower edges of the box are 75%，the upper and lower limits are 90%，and the outer shape is the kernel density estimation.圖1 FPKM密度分布Figure 1 Density distribution of FPKM

2.4 差異表達基因功能富集分析

2.4.1 GO富集分析 對差異表達基因進行GO富集，結果顯示，與細胞組分(CC)相關的GO條目有491個，占總量的12.68%；與分子功能(MF)相關的GO條目有1 062個，占總量的27.42%；與生物過程(BP)相關的GO條目有2 320個，占總量的59.90%.對GO富集分析結果，每種分類挑選P-value最小，即富集最顯著的前10個GO條目進行展示，結果見圖4.CC分類中，P-value值最小的基因的ID號為GO:0005887，與質膜組成成分相關；MF分類中，P-value值最小的基因的ID號為GO:0004046，與氨酰基轉移酶活性相關；BP分類中，P-value值最小的基因的ID號為GO:0044281，與小分子代謝過程相關.

x表示差異倍數；圖中兩條豎虛線為2倍表達差異閾值；橫虛線為P值的閾值(0.05).紅點表示該組上調基因，藍點表示該組下調基因，灰點表示非顯著差異表達基因.x means fold change；The two vertical dashed lines in the figure are 2 times of the expression difference threshold；the horizontal dashed line is the threshold value of P-value(0.05).The red，blue and gray dots represent the up-regulated，down-regulated and non-significantly differentially expressed genes，respectively.圖2 差異表達基因的火山圖Figure 2 Volcano map of differentially expressed genes

2.4.2 KEGG富集分析 參考KEGG通路系統的基因組信息數據庫，對檢測到的在這兩個玉米花粉材料中呈現差異表達的基因的代謝通路作富集分析.富集到的差異基因數目為349個，共涉及到102種通路，這些小通路可進一步歸為6個大通路，包括細胞過程3個、環境信息處理4個、遺傳信息處理20個、人類疾病1個、新陳代謝72個和有機系統2個.上調和下調表達基因分別有143和206個(注：一個基因可能同時富集于不同通路).挑選出的P-value最小即富集最顯著的前30個Pathway展示結果見圖5，其中1個與環境信息處理有關，4個與遺傳信息處理有關，25個與新陳代謝有關.

最外圈是染色體條帶.紅色和綠色分別為上調和下調基因的log2x(Fold Change) 值的柱狀圖，灰色為無差異表達基因的log2x(Fold Change)值的散點圖.The outermost circle is the chromosomal banding.The red and green are the histogram of the log2x(Fold Change) value of up-regulated and down-regulated genes respectively，and the gray is the scatter graph of the log2x(Fold Change) value of the undifferentially expressed genes.圖3 基因組圈圖Figure 3 Genome circos

3 討論與結論

溫度是影響植物生長非常重要的因素[12-13]，高溫脅迫是最常見的非生物脅迫之一，它能顯著抑制植物生長發育[14].高溫脅迫一般會降低水分含量，影響植物光合作用和呼吸作用[15].研究發現，高溫脅迫已導致全球植物的產量下降和品質降低[16-17].

1：質膜組成成分;2：DNA指導的RNA聚合酶II，全酶；3：胞質；4：核小體；5：胞質大核糖體亞基；6：DNA包裝復合物；7：核DNA指導RNA聚合酶復合物；8：大核糖體亞基；9：MLL1/2復合物；10：MLL1復合物；11：氨基酰化酶活性；12：蛋白異源二聚化活性；13：基本轉錄機制結合；14：基本RNA聚合酶II轉錄機制結合；15：肽基轉移酶活性；16：谷胱甘肽水解酶活性；17：輔因子結合；18：核苷跨膜轉運體活性；19：磷酸吡哆醛結合；20：維生素B6結合；21：小分子代謝過程；22：赤霉素生物合成過程；23：應答低氧含量；24：應答氧水平；25：糖脂分解代謝過程；26：鞘糖脂分解代謝過程；27：萜類生物合成過程；28：磷代謝過程的正向調節過程；29：磷酸鹽代謝過程的正向調節過程；30：有機酸代謝過程.1：Integral component of plasma membrane；2：DNA-directed RNA polymerase II，homoenzyme；3：Cytosol；4：Nucleosome；5：Cytosolic large ribosomal subunit；6：DNA packaging complex；7：Nuclear DNA-directed RNA polymerase complex；8：Large ribosomal subunit；9：MLL1/2 complex；10：MLL1 complex；11：Aminoacylase activity；12：Protein heterodimerization activity；13：Basal transcription machinery binding；14：Basal RNA polymerase II transcription machinery binding；15：Peptidyltrasferase activity；16：Glutathione hydrolase activity；17：Cofactor binding；18：Nucleoside transmembrane transporter activity；19：Pyridoxal phosphate binding；20：Vitamin B6 binding；21：Small molecule metabolic process；22：Gibberellin biosynthetic process；23：Response to decreased oxygen levels；24：Response to oxygen levels；25：Glycolipid catabolic process；26：Glycosphingolipid catabolic process；27：Terpenoid biosynthetic process；28：Positive regulation of phosphorus metabolic process；29：Positive regulation of phosphate metabolic process；30：Organic acid metabolic process.圖4 差異表達基因的GO注釋和分類Figure 4 GO annotation and classification of differentially expressed genes

玉米是世界上最重要的農作物之一，高溫脅迫同樣會顯著降低其產量和品質[18].因此，研究玉米對高溫脅迫響應的分子機制及選育耐高溫脅迫的玉米品種具有重要意義.

雖然已有關于玉米耐高溫脅迫的基因克隆或表達方面的研究，但是以前關于基因克隆的研究主要限于單個或少數幾個基因的克隆表達分析，如Wang等[19]從玉米中克隆了轉錄因子ZmWRKY106基因，將其轉入擬蘭芥，過量表達的結果顯示該基因會明顯提高對干旱和高溫脅迫的抗性；Ma等[20]研究發現ZmbZIP4基因的過表達可以調節ABA的生物合成及根的發育，從而提高玉米包括高溫脅迫等在內的多種抗性；Li等[21]也發現玉米ZmbZIP60基因的表達會應答于高溫脅迫.由于玉米對高溫脅迫的抗逆性屬多基因控制的數量性狀，因而對單個基因的克隆表達研究還不能完全滿足育種需要，有必要在基于轉錄組測序的在全基因組水平下研究耐高溫脅迫基因.

基于轉錄組測序對玉米耐高溫脅迫的基因差異表達以前也有類似研究，且發現了大量與高溫脅迫相關的差異表達基因，但以前的研究取材部位主要限于玉米葉片[1，18，22-23].然而在包括中國在內的許多國家，在玉米散粉期常遭遇高溫天氣，高溫脅迫導致花粉活力下降，從而影響玉米的授粉、結實，使其大幅減產，因此，選擇以花粉為材料，基于轉錄組高通量測序，研究不同玉米材料間的基因差異表達對于認識玉米耐高溫脅迫的分子遺傳機制具有重要意義.

本試驗在研究方法、材料選擇或取材部位等方面，不同于以前的研究[19-23]，本試驗是以中地88(耐高溫脅迫)和先玉335(高溫敏感)兩個玉米品種的花粉為研究對象，研究其基因表達情況，共發現差異表達基因1 822個，相對于對照(先玉335)而言，中地88上調表達基因有857個，下調表達基因965個.另外還發現有22 731個基因在兩個材料中表達差異不顯著.對差異表達基因進行GO功能注釋，結果顯示大部分基因都與細胞組分、分子功能或生物過程功能相關.KEGG功能富集結果顯示，這些差異表達基因可以歸為6個Level 1分類(包括102個Levle 2通路).綜合分析這些差異表達基因，發現有7個基因直接應答于熱激效應，其中2個上調表達，5個下調表達；有3個基因與熱激蛋白結合有關，其中上調表達基因1個，下調表達基因2個.在這12個與玉米耐高溫脅迫相關的基因中，多數為在玉米中新發現的與耐高溫脅迫相關的基因.依據基因的ID號和玉米參考基因組序列信息，參考以前的研究[24-26]，進一步對其基因全長序列、cDNA序列及表達蛋白的生物學特性等作進一步研究，關于這些基因的分子結構及遺傳功能研究正在進行之中.本實驗的研究結果為進一步闡明玉米耐高溫脅迫的分子遺傳機制奠定了一些理論基礎，為在植物中克隆應用這些耐高溫脅迫基因提供了參考.

1：植物MAPK信號通路；2：核糖體；3：堿基切除修復；4：硫傳遞系統；5：RNA運輸；6：其他多糖降解；7：糖酵解/糖異生；8：α-亞麻酸代謝；9：淀粉和蔗糖代謝；10：苯并噁唑嗪酮類化合物生物合成；11：鞘脂代謝；12：丙酮酸代謝；13：丁酸代謝；14：脂肪酸降解；15：磷酸戊糖途徑；16：氰氨基酸代謝；17：抗壞血酸和醛糖酸代謝；18.苯丙氨酸代謝；19：咖啡因代謝；20：纈氨酸、亮氨酸和異亮氨酸降解；21：組氨酸代謝；22：苯丙氨酸、酪氨酸、色氨酸生物合成；23：果糖和甘露糖代謝；24：亞油酸代謝；25：牛磺酸與次牛磺酸代謝；26：煙酸和煙酰胺代謝；27：半乳糖代謝；28：硫胺代謝；29：有機含硒化合物代謝；30：泛醌和其他萜醌生物合成1：MAPK signaling pathway - plant；2：Ribosome；3：Base excision repair；4：Sulfur relay system；5：RNA transport；6：Other glycan degradation；7：Glycolysis / Gluconeogenesis；8：alpha-Linolenic acid metabolism；9：Starch and sucrose metabolism；10：Benzoxazinoid biosynthesis；11：Sphingolipid metabolism；12：Pyruvate metabolism；13：Butanoate metabolism；14：Fatty acid degradation；15：Pentose phosphate pathway；16：Cyanoamino acid metabolism；17：Ascorbate and aldarate metabolism；18：Phenylalanine metabolism；19：Caffeine metabolism；20：Valine，leucine and isoleucine degradation；21：Histidine metabolism；22：Phenylalanine，tyrosine and tryptophan biosynthesis；23：Fructose and mannose metabolism；24：Linoleic acid metabolism；25：Taurine and hypotaurine metabolism；26：Nicotinate and nicotinamide metabolism；27：Galactose metabolism；28：Thiamine metabolism；29：Selenocompound metabolism；30：Ubiquinone and other terpenoid-quinone biosynthesis圖5 KEGG通路富集結果柱狀圖Figure 5 Histogram of enrichment results for KEGG pathway