唐 軍， 丁西朋， 陳志堅， 馬向麗， 畢玉芬， 郭鳳根*

(1. 云南農(nóng)業(yè)大學動物科學技術學院， 云南 昆明 650201； 2. 中國熱帶農(nóng)業(yè)科學院熱帶作物品種資源研究所/農(nóng)業(yè)部華南作物基因資源與種質創(chuàng)制重點實驗室， 海南 海口 571101)

木豆為世界第六大食用豆類，多年生常綠小灌木，株高1.5～3.0 m。木豆的用途非常廣泛，嫩莖葉用作動物飼料或藥用，還可用作水土保持和覆蓋作物，主要分布在世界熱帶和亞熱帶地區(qū)[1-2]。在我國主要分布于長江以南的省(區(qū))，包括海南、廣西、廣東、云南等地[3]。低溫脅迫是作物生長過程中遭受的重要環(huán)境脅迫因素之一，可以引起在細胞膜系統(tǒng)、酶活性、細胞質失水等，進而細胞代謝失衡或衰亡等[4-5]。南方地區(qū)時常受到低溫的影響，導致木豆籽實產(chǎn)量變低，甚至死亡，對木豆的生產(chǎn)受到了限制。深入研究木豆的耐寒機制，對后期篩選優(yōu)異的耐寒種質資源和新品種的培育，解決南方木豆生產(chǎn)受低溫影響具有重要的意義。

轉錄組是生物體特定組織或特定細胞在某一階段產(chǎn)生的所有轉錄的集合，可以獲取細胞內所有基因轉錄后的信息，可以對基因結構及表達進行深入和更詳細的分析[6-7]。轉錄組在植物低溫脅迫條件上的研究的應用較多，已成為植物抗逆研究等重要的手段。筆者在查閱利用轉錄組測序技術揭示木豆耐低溫機制方面研究文獻，發(fā)現(xiàn)國內外少有報道。為研究木豆響應低溫脅迫的機制，在前期試驗獲得1份耐低溫和1份敏感試驗材料基礎上，進行低溫脅迫處理，利用轉錄組進行了研究，以期從分子層面上解析木豆耐低溫的機制，為木豆的耐寒育種提供基礎理論資料。

1 材料與方法

1.1 試驗材料與處理條件

試驗選用耐低溫緬甸木豆(MD)和低溫敏感型‘瓊中木豆’(QZ)材料各1份，由中國熱帶農(nóng)業(yè)科學院熱帶作物品種資源研究所提供。將‘緬甸木豆’、瓊中木豆種子直播于裝有本地磚紅壤的育苗袋(10 cm×10 cm)中育苗，每袋5粒。放置溫室大棚中30 d左右，待幼苗長至5～7葉。將2份材料置于低溫庫(4℃，光照10 h，濕度75%～80%)中處理96 h。將MD，QZ木豆對照和處理(MDCK/MDT，QZCK/QZT)的葉片分別進行混合取樣。每個處理設3個生物學重復，并進行編號，共計12份樣品，送華大基因公司DNBSEQ平臺進行轉錄組測序。

1.2 RNA提取、cDNA文庫構建

總RNA采用常規(guī)CTAB法進行提取。cDNA文庫構建流程：mRNA分離將mRNA片段化合成二鏈cDNA末端修復和接頭連接PCR反應及產(chǎn)物回收PCR產(chǎn)物環(huán)化文庫檢測PCR產(chǎn)物環(huán)化上機測序。

1.3 測序數(shù)據(jù)的過濾、組裝和基因表達量計算

將檢測數(shù)據(jù)使用Soapnuke(v1.4.0)進行過濾，得到Clean Data，利用Hisat軟件將clean reads比對木豆參考基因組(版本：GCF_000340665.1_C.cajan_V1.0_NCBI[8])進行序列比對，對得到的Mapped Data進行插入片段長度檢驗、隨機性檢驗等文庫質量評估。使用Bowtie2(v2.2.5) 軟件將clean reads比對到參考基因序列，之后再使用RSEM(v1.2.8)計算基因和轉錄本的表達水平。轉錄因子使用getorf(EMBOSS:6.5.7.0)檢測Unigene的ORF，并鑒定。

1.4 差異表達基因(Differentially expressed genes，DEGs)、GO與KEGG注釋和富集

DEGs篩選通過DEseq2方法基于負二項分布原理方法進行DEGs檢測，參數(shù)：Qvalue≤0.05，并進行GO功能和KEGG通路注釋。同時使用R軟件中的phyper函數(shù)進行富集分析，Qvalue≤0.05的功能為顯著富集。

1.5 目標基因qRT-PCR驗證

為深入了解測序數(shù)據(jù)的正確性，本研究選擇了10個可能在木豆響應低溫脅迫過程中發(fā)揮作用的差異基因進行qRT-PCR驗證。使用Qusant StudioTMReal-Time PCR(Thermo Fisher，美國)系統(tǒng)和SYBR Green(Vazyme，中國)定量試驗盒進行實時熒光定量PCR反應。選擇木豆IF4α基因為內參基因[9]，每個樣品3次重復，以2-ΔΔCt法計算差異基因的表達水平。試驗引物如表1所示。

表1 實時熒光定量PCR基因和引物Table 1 Genes and primers used for qRT-PCR

1.6 數(shù)據(jù)分析和圖處理

EXCEL2019、R語言(version 4.1.0)和華大數(shù)據(jù)庫分析軟件處理相關數(shù)據(jù)和繪圖。

2 結果與分析

2.1 轉錄組測序和組裝

2.1.1轉錄組數(shù)據(jù)質控分析 對MDCK、MDT、QZCK和QZT的木豆幼苗葉片進行RNA-seq測序，測序結果顯示如表2所示，12個樣品共獲得76.75 Gb Clean Data，各樣本平均產(chǎn)出超過6.40 Gb的高質量數(shù)據(jù)，Q20堿基百分比均在95.78%及以上，Q30堿基百分比均在90.16%及以上，樣本與參考基因比對比例89.22%～91.97%，說明測序過程中建庫質量較好，如表2所示。

表2 低溫脅迫下木豆葉片樣本轉錄組數(shù)據(jù)統(tǒng)計Table 2 Statistics of transcriptome sequencing of pigeonpea leaf samples under cold stress

2.1.2轉錄組測序結果 樣品比對基因組的平均比對率為90.42%，比對基因集的平均比對率為82.42%；共檢測到表達的基因數(shù)為28 796個，其中已知的基因為26 452個，預測的新基因為2 344個；共檢測出21 760個新轉錄本，其中15 326個屬于已知蛋白編碼基因的新的可變剪接亞型，2 474個屬于新的蛋白編碼基因的轉錄本，剩下的3 960個屬于長鏈非編碼RNA。

2.2 轉錄因子鑒定

本研究注釋的轉錄因子1 916個，共計58類基因家族，其中MYB家族轉錄因子占14.43%，AP2-EREBP家族轉錄因子占9.48%，bHLH家族轉錄因子占8.81%，WRKY家族轉錄因子占5.21%，說明這四類轉因子在木豆響應低溫脅迫中占有主導地位，前20類轉錄因子如表3示。

表3 低溫脅迫下木豆的主要轉錄因子鑒定Table 3 Identification of transcript factors in pigeonpea under cold stress

2.3 差異基因的表達分析

2.3.1樣品間差異基因表達量分析 相關性熱圖顯示MD、QZ兩份材料重復間皮爾遜相關系數(shù)較高，而參試材料間相關系數(shù)較小，說明樣品處理較好。

圖1 12個木豆樣本間的FPKM皮爾遜系數(shù)相關性熱圖Fig.1 The heat map of Pearson correlation of FPKM in 12 pigeonpea samples

2.3.2組間差異表達基因分析 MDT-MDCK間差異表達基因數(shù)達16 407個(上調8 538，下調7 869)；QZT-QZCK間差異表達基因有16 021個(上調8 471，下調7 550)。MDT-QZT和MDCK-QZCK差異基因分別為6 027個(上調3 440，下調2 587)和5 994個(上調3 456，下調2 538)。耐低溫材料MD較低溫敏感材料QZ的差異基因多15.57%，說明耐低溫材料的差異基因表達更為豐富，如圖2所示。

圖2 MDT-MDCK,QZT-QZCK,MDT-QZT和MDCK-QZCK組間差異表達基因Fig.2 Number of DEGs from MDT-MDCK,QZT-QZCK,MDT-QZT and MDCK-QZCK

MDT-MDCK,QZT-QZCK,MDT-QZT和MDCK-QZCK差異基因為20 034個。MDT-MDCK差異表達基因占81.90%，QZT-QZCK差異表達基因占79.97%，MDT-QZT差異基因占30.08%，MDCK-QZCK差異基因占29.92%，同樣說明耐低溫材料MD中顯著差異表達基因較多，如圖3示。同時，分析發(fā)現(xiàn)4個組間同步上調表達的差異基因為4 121個，同步下調表達的差異基因4 241個。

圖3 MDCK-MDT,QZCK-QZT，QZT-MDT和MDCK-QZCK差異表達基因維恩圖Fig.3 Venn diagram of DEGs from MDCK-MDT,QZCK-QZT，QZT-MDT and MDCK-QZCK

2.4 GO注釋和富集分析

GO功能分析將表達基因分為三個大類，即生物學過程、分子功能和細胞組成。通過GO功能注釋，MDCK-MDT組和QZCK-QZT組將差異表達基因歸于生物學過程(13個亞類)、分子功能(14個亞類)和細胞組成(21個亞類)共計48個亞類，且均在分子功能過程中“細胞”、“膜”、細胞組分部分中“結合”和“催化活性”、生物學過程中“細胞過程”、“代謝過程”等功能中注釋基因相對較多，能被注釋功能的差異基因分別為13 524個和13 249個，MD較QZ多出2.08%。MDCK-MDT組在分子功能中“結合”、“催化活性”中差異基因有6 659個(上調3 669，下調2 990)和6 651個(上調3 306，下調3 345)。細胞組分中“細胞”、“膜”中差異基因分別為6 294個(上調3 357，下調2 937)和4 982個(上調2 281，下調2 701)。生物學過程中“細胞過程”和“代謝過程”注釋的基因相對較多，分別有5 079個(上調2 753，下調2 326)和4 334個(上調2 323，下調2 011)。前20條GO注釋如圖4 A所示。

QZCK-QZT組在生物學過程中，分子功能中“結合”、“催化活性”中有6 539個(上調3 608，下調2 931)和6 494個(上調3 234，下調3 260)基因。細胞組分中“細胞”、“膜”注釋的差異基因分別為6 176個(上調3 358，下調2 818)和4 877個(上調2 260，下調2 617)。生物學過程“細胞過程”、“代謝過程”富集較多，分別為5 001個(上調2 268，下調2 733)和4 309個(上調2 338，下調1 971)。前20條GO注釋，如圖4 B所示。

圖4 低溫脅迫下不同組間差異表達基因GO功能分析(A. MDT-MDCK;B. QZT-QZCK)Fig.4 GO classification analysis of DEGs under cold stress(A. MDT-MDCK;B. QZT-QZCK)注：Ⅰ. 分子功能(1～5);Ⅱ. 細胞組成(6～12);Ⅲ. 生物學過程(13～20)。1. 結合；2. 催化活性；3. 轉運活性；4. 轉錄調節(jié)活性；5.結構分子活性；6. 細胞；7膜;8細胞器；9. 膜部件；10. 細胞器部分；11. 含蛋白質復合物；12. 膜腔；13. 細胞過程；14. 代謝過程；15.生物調節(jié)；16. 應激響應；17. 細胞成分組織或生物發(fā)生；18. 定位；19. 信號；20. 發(fā)展過程Note：Ⅰ. Molecular function(1～5);Ⅱ. Cellular component(6～12);Ⅲ. Biological processs(13～20). 1. Binding；2. Catalytic activity；3. Transporter activity；4. Transcription regulator activity；5. Structural molecule activity；6. Cell；7. Membrane;8. Organelle；9. Membrane part；10. Organelle part；11. Protein-containing complex；12. Membrane-enclosed lumen；13. Cellular process；14. Metabolic process；15. Biological regulation；16. Response to stimulus；17. Cellular component organization or biogenesis；18. Localization；19. Signaling；20. Developmental process

GO富集分析發(fā)現(xiàn)，低溫脅迫誘導MDT-MDCK上調基因顯著富集GO分類條目194條，前2條為在細胞組成中“細胞”(GO:0005622)和生物學過程中“細胞過程”(GO:0043604)，下調基因顯著富集GO分類條目132條，前2條為細胞組成中“膜”(GO:0016020)和生物學過程中“細胞過程”(GO:0006464)。低溫脅迫誘導QZT-QZCK上調基因最顯著富集在生物過程中“細胞過程”(GO：0043604；GO:0006412)，下調基因顯著富集在細胞組成中“膜”(GO:0016020)和生物學過程中“細胞過程”(GO:0006464)。

2.5 KEGG途徑分類和富集分析

2.5.1KEGG途徑分類 MD和QZ兩份材料KEGG注釋到分別有12 046個和12 967個DEGs，數(shù)量由高到低分為代謝過程、遺傳信息處理、環(huán)境信息處理、細胞過程和有機系統(tǒng)(含19個代謝通路)。MDT-MDCK的KEGG注釋到DEGs分別為代謝過程中的“全球和概覽圖”DEGs為3 008個(上調1 573，下調1 435)，遺傳信息處理中“翻譯”1 064個(上調640，下調424)、有機系統(tǒng)中的“環(huán)境適應”934個(上調503，下調431)、環(huán)境信息處理中的“信號轉導”643個(上調254，下調389)、細胞過程中“運輸和分解”592個(上調254，下調338)。如圖5A所示。

QZT-QZCK的DEGs分別為代謝過程中的“全球和概覽圖”2 978個(上調1 586，下調1 392)、遺傳信息處理中“翻譯”1 086個(上調665，下調421)、有機系統(tǒng)的“環(huán)境適應”939個(上調486，下調453)，環(huán)境信息處理中的“信號轉導”651個(上調402，下調249)、細胞過程中“運輸和分解”573個(上調320，下調253)。如圖5B所示。

圖5 低溫脅迫下不同組間差異表達基因KEGG途徑(A.MDT-MDCK;B.QZT-QZCK)Fig.5 KEGG pathway of different expression gene under cold stress (A. MDT-MDCK;B.QZT-QZCK)注：Ⅰ.有機系統(tǒng)(1);Ⅱ.代謝(2～12);Ⅲ.遺傳信息處理(13～16);Ⅳ.環(huán)境信息處理(17,18);Ⅴ.細胞過程(19)。1.環(huán)境適應;2.全球和概覽圖;3.碳水化合物代謝;4.脂質代謝;5.其他次生代謝產(chǎn)物的生物合成;6.氨基酸代謝;7.能量代謝;8.輔助因子和維生素的代謝;9.其他氨基酸代謝;10.聚糖生物合成與代謝;11.萜類化合物和聚酮類化合物的代謝;12.核苷酸代謝;13.翻譯;14.折疊、分類和降解;15.轉錄;16.復制和修復;17.信號轉導;18.膜運輸;19.運輸和分解Note：Ⅰ. Organismal Systems (1);Ⅱ. Metabolism (2～12);Ⅲ. Genetic Information Processing (13～16);Ⅳ. Environmental Information Processing (17,18);Ⅴ. Cellular Processes (19). 1. Environmental adaptation;2. Global and overview maps;3. Carbohydrate metabolism;4. Lipid metabolism;5. Biosynthesis of other secondary metabolites;6. Amino acid metabolism;7. Energy metabolism;8. Metabolism of cofactors and vitamins;9. Metabolism of other amino acids;10. Glycan biosynthesis and metabolism;11. Metabolism of terpenoids and polyketides;12. Nucleotide metabolism;13. Translation;14. Folding,sorting and degradation;15. Transcription;16. Replication and repair;17. Signal transduction;18. Membrane transport;19. Transport and catabolism

2.5.2KEGG通路富集分析 低溫脅迫后，經(jīng)KEGG通路分析，差異基因主要富集的前10條KEGG途徑如下表所示。MDT-MDCK組顯著富集在5個代謝通路，即核糖體、植物激素信號轉導通路、植物晝夜節(jié)律通路、光合作用天線蛋白通路、卟啉與葉綠素代謝通路，其中植物激素信號傳導通路最多，差異基因富集數(shù)量達543個，其次為核糖體通路為341個基因(表4)。QZT-QZCK組差異基因顯著富集在3個代謝途徑，即核糖體通路、磷酸肌醇代謝和磷脂酰肌醇信號系統(tǒng)，通路富集基因數(shù)分別為353，119和113個；富集差異基因最多的為植物激素信號傳導通路，達513個，其次為植物促分裂素原活化蛋白激酶(MAPK)信號通路，差異基因達387個，差異不顯著(表5)。耐低溫材料MD中顯著富集的KEGG通路較低溫敏感材料QZ中多，兩者均共同富集的通路為核糖體通路。植物激素信號轉導通路、植物晝夜節(jié)律通路、光合作用天線蛋白通路、卟啉與葉綠素代謝通路只在MD材料中顯著富集，說明這些通路可能參與耐寒材料響應低溫脅迫的表達。

表4 MDT-MDCK差異表達基因KEGG途徑Table 4 KEGG pathway of DEGs in MDT-MDCK

表5 QZT-QZCK差異表達基因KEGG途徑Table 5 KEGG pathway of DEGs in QZT-QZCK

2.6 耐低溫脅迫響應的差異基因分析

耐寒基因在耐寒材料中顯著表達，敏感材料中表達差異不顯著的特點。由圖3知，MDT-MDCK、MDT-QZT共同表達有基因528個，設置差異基因|log2Fold Change|≥1，顯著富集Q≤0.05，QZCK-QZT和MDCK-QZCK的Q≥0.05為篩選條件，分析發(fā)現(xiàn)4個組間DEGs為150個，其中上調DEGs為98個，下調DEGs為52個。上調DEGs功能顯著性富集在分子功能中“催化活性”中α-1,4-葡萄糖苷酶(GO:0004558)和麥芽糖α-葡萄糖苷酶(GO:0032450)。KEGG通路分析發(fā)差異基因在代謝中的“全球和概覽圖”、“碳水化合物代謝”和“脂質”最多。KEGG富集分析發(fā)現(xiàn)代謝中“脂質代謝”(ko04016)和“次級代謝物合成”(ko00940)。下調DEGs功能顯著富集在分子功能中“催化活性”(GO:0004558)和“結合”(GO:0005488)。KEGG通路分析發(fā)差異基因在代謝中的“全球和概覽圖”、“碳水化合物代謝”最多，顯著富集代謝中“碳水化合物代謝”(ko00500)、“聚糖生物合成與代謝”(ko00511)和“脂質代謝”(ko04016)。結果說明這些基因可能參與在耐寒材料中響應低溫脅迫表達，在木豆耐寒特性中發(fā)揮著重要作用。

2.7 其他與低溫脅迫相關的差異基因及通路

基于相關文獻，通過MDT-MDCK、QZT-QZCK組間差異表達基因和KEGG富集通路比對，發(fā)現(xiàn)了37個基因和14個KEGG通路在木豆響應中差異較大，這些差異基因和通路也可能與木豆耐低溫的特性相關，如表6所示。

表6 木豆響應低溫脅迫差異基因Table 6 DEGs in pigeonpea responses to cold tolerance

2.6 實時熒光定量PCR驗證

選取10個差異基因進行驗證，在MD和QZ中的差異表達與轉錄組基因定量表達趨勢相同，證明本次測序可信度高，為后續(xù)開展耐寒基因深入挖掘提供詳實的信息。

圖6 差異表達基因的qRT-PCR驗證Fig.6 Validation of DEGs by qRT-PCR

3 討論

低溫對植物生理和代謝過程受溫度差異、脅迫時間、基因型等影響存在差異，低溫導致光合同化產(chǎn)物積累減少，同時降低了能量傳輸，降低生化反應酶的活性，對植物的代謝平衡有著嚴重的影響10]。

受低溫脅迫后，木豆耐低溫材料的差異基因多于低溫敏感材料。耐低溫脅迫下青海早熟禾[11]、茄子[12]、苦瓜[13]、萱草[14]、花生[15]和水稻[16]等牧草和作物中轉錄組研究結果也發(fā)現(xiàn)耐低溫材料中的差異基因要多于敏感材料。研究結果一定程度上說明耐寒型相較于低溫敏感型木豆在低溫環(huán)境時，有更多的特異基因表達來抵御低溫逆境的影響。

轉錄調控是植物對逆境脅迫產(chǎn)生應答的關鍵步驟，轉錄因子在植物對逆境脅迫的應答過程中發(fā)揮重要的作用[17-18]。有研究發(fā)現(xiàn)MYB轉錄因子參與植物對低溫脅迫的響應[19]。枳MYB轉錄因子家族PtsrMYB，該基因的轉錄水平被非生物調控如脫水、鹽、冷和ABA處理表達上調[20]。番茄MYB15在低溫抗性的CBF途徑中起正調控作用[21]。甘薯中有36個R2R3-MYB轉錄因子基因可能在響應生物和非生物脅迫上發(fā)揮重要功能[22]。水稻MYB轉錄因子在的冷適應中起著關鍵調控作用[23]?；ㄉ劝l(fā)期響應低溫脅迫轉錄因子基因bHLH、MYB和EREB的表達量明顯升高，F(xiàn)ER在脅迫12 h后也有明顯升高[15]。對水稻耐寒材料的轉錄組分析發(fā)現(xiàn)有1654個基因特異性地改變了表達冷應激，基因主要屬于轉錄因子家族AP2-EREBP、NAC和WRKY[24]。從小麥基因組中克隆得到TaWRKY10基因，在聚乙二醇、NaCl、低溫和H2O2處理后上調[25]。金銀花12個WRKY基因與已知抗寒基因有同源性，參與低溫逆境響應[26]。Singh等[27]研究發(fā)現(xiàn)WRKY轉錄因子參與木豆逆境下的表達。黃花苜蓿響應低溫脅迫的轉錄因子中WRKY、ERF、MYB、bHLH和NAC等5個轉錄因子家族包含的差異表達基因數(shù)目較多[28]。低溫脅迫紫花苜蓿中差異表達基因中發(fā)現(xiàn)了43個轉錄因子家族，共251個基因被顯著地誘導表達，且轉錄因子家族基因受低溫脅迫誘導表達不同[29]。花煙草低溫脅迫下有118個轉錄因子的表達發(fā)生了顯著改變，包括NAC、ERF、MYB、WRKY等[30]。在本文研究中檢測到木豆中MYB、AP2-EREBP、bHLH和WRKY等轉錄因子，與前人在研究中植物響應低溫脅迫中的發(fā)現(xiàn)類似。

低溫脅迫后，GO分析表明DEGs主要注釋在分子功能中“結合”和“催化活性”、細胞組成次級類別中“細胞”及“膜”系統(tǒng)、生物學過程中“細胞過程”、“代謝過程”的差異基因較多，說明對于木豆受低溫脅迫后，木豆分子結合功能、催化活性、細胞結構及膜系統(tǒng)受損，并影響細胞過程及代謝過程，進而導致差異表達基因數(shù)量增加，形成一個抵御逆境復雜過程。KEGG分析發(fā)現(xiàn)在全球和概覽圖、翻譯、環(huán)境適應、信號轉導、細胞過程中運輸和分解代謝通路富集到的DEGs數(shù)目較多。進一步富集分析發(fā)現(xiàn)耐寒MD材料中差異基因顯著富集在植物激素信號轉導通路、植物晝夜節(jié)律通路、光合作用天線蛋白通路、卟啉與葉綠素代謝通路等，說明這幾類通路在耐寒木豆材料中發(fā)揮重要作用，可以進一步進行深入研究。相關類似報道如野生早熟禾、苦瓜、花生低溫脅迫轉錄組分析結果相似[11,13,15]。低溫能抑制了植物光合作用，且破壞光能吸收、轉運的進行[31]。茶樹低溫響應中，差異表達基因大量富集在葉綠體和光合作用上[32]。油菜轉錄組測序發(fā)現(xiàn)抗寒品系響應冷凍脅迫的表達基因富集在“植物激素信號轉導”、“碳水化合物代謝”、“氨基酸代謝”、“脂質代謝”等代謝途徑中的差異基因在冷凍脅迫下被誘導表達[33]。川百合寒冷應答的分子機制經(jīng)富集分析后與光合作用、光合作用-天線蛋白等代謝相關[34]。東鄉(xiāng)野生稻的光合作用代謝參與響應低溫脅迫過程[35]。百香果低溫脅迫轉錄組植物激素信號轉導途徑參與百香果響應低溫脅迫[36]。蘋果花芽早期響應低溫信號的基因表達情況及碳水化合物有關的代謝生物學過程，DEGs主要參與植物激素信號轉導等代謝途徑[37]。耐寒花生種質其中植物激素信號轉導、淀粉和蔗糖代謝、植物晝夜節(jié)律和α-亞麻酸代謝途徑相關[38]。低溫脅迫對水稻的影響主要體現(xiàn)在生物學過程中，在光合作用、次生代謝的生物合成和苯丙氨酸代謝均有明顯作用[39]。低溫脅迫冰菜植物激素信號轉導、淀粉和蔗糖代謝等途徑負向影響[40]。這些通路有可能是木豆耐低溫材料響應低溫脅迫特異表達通路，與耐寒特性相關。MD和QZ兩份材料中GO分析DEGs參與細胞組成中“膜”和生物學過程中“細胞過程”。KEGG分析MD材料富集的通路不同，MD較QZ多，MD主要富集在植物激素信號轉導通路、核糖體等通路中。在MD材料中，發(fā)現(xiàn)150個顯著差異DEGs，參與碳水化合物、脂質代謝途徑，與木豆耐寒特性密切相關。木豆響應低溫脅迫的過程主要通過促調節(jié)糖分解酶類活性來促進糖的分解、脂質代謝增加細胞結構的穩(wěn)定性，緩解傷害。

植物細胞膜感受到低溫后，啟動信號分子傳遞，然后調控相關基因表達，最終緩解低溫危害[14]。另有研究報道植物激素信號傳導、酪氨酸代謝、脯氨酸代謝、MAPK信號通路、Ca2+信號通路等參與了植物逆境信號傳導[11,41-45]。本試驗中耐寒型MD木豆材料也檢測到上述代謝通路差異基因表達存在表達，說明這些信號通路均響應耐寒脅迫。研究發(fā)現(xiàn)熱休克蛋白(HSPs)參與木豆對低溫脅迫的響應，這與萱草HSP蛋白基因冷敏感型和耐冷型萱草中呈現(xiàn)差異表達[14]、報春花熱應激HSP蛋白基因冷脅迫條件下也有表達[46]。信號通路、能量代謝通路及敏感基因均參與響應，說明木豆對低溫脅迫是一個復雜的過程，耐寒的表型是綜合系統(tǒng)響應的過程。

4 結論

本研究以耐低溫和低溫敏感木豆為研究材料，通過RNA-seq獲得了豐富的轉錄本信息，對差異表達基因的基因功能進行注釋和分析，MD和QZ兩份材料中GO分析發(fā)現(xiàn)DEGs大部分注釋在細胞組成中“膜”和生物學過程中“細胞過程”，KEGG分析MD材料富集的通路不同，MD較QZ多，MD主要富集在植物激素信號轉導通路、核糖體等通路中。木豆響應低溫脅迫的轉錄因子主要為主要包括MYB，AP2-EREBP等。兩份材料比較發(fā)現(xiàn)在MD中特異的150個DEGs，主要為糖分解酶類，參與碳水化合物、脂質代謝途徑，可能與耐低溫的特性密切相關。同時篩選到37個與耐寒相關重要差異表達基因和14條KEGG通路。本文從分子水平闡述了木豆低溫脅迫下的響應機制，為后期木豆耐寒育種提供了參考資料。