木豆響應(yīng)低溫脅迫差異表達基因分析

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

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

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

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

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)的葉片分別進行混合取樣。每個處理設(shè)3個生物學重復，并進行編號，共計12份樣品，送華大基因公司DNBSEQ平臺進行轉(zhuǎn)錄組測序。

1.2 RNA提取、cDNA文庫構(gòu)建

總RNA采用常規(guī)CTAB法進行提取。cDNA文庫構(gòu)建流程：mRNA分離將mRNA片段化合成二鏈cDNA末端修復和接頭連接PCR反應(yīng)及產(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進行插入片段長度檢驗、隨機性檢驗等文庫質(zhì)量評估。使用Bowtie2(v2.2.5) 軟件將clean reads比對到參考基因序列，之后再使用RSEM(v1.2.8)計算基因和轉(zhuǎn)錄本的表達水平。轉(zhuǎn)錄因子使用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個可能在木豆響應(yīng)低溫脅迫過程中發(fā)揮作用的差異基因進行qRT-PCR驗證。使用Qusant StudioTMReal-Time PCR(Thermo Fisher，美國)系統(tǒng)和SYBR Green(Vazyme，中國)定量試驗盒進行實時熒光定量PCR反應(yīng)。選擇木豆IF4α基因為內(nèi)參基因[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ù)庫分析軟件處理相關(guān)數(shù)據(jù)和繪圖。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測序和組裝

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

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

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

2.2 轉(zhuǎn)錄因子鑒定

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

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

2.3 差異基因的表達分析

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

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

2.3.2組間差異表達基因分析 MDT-MDCK間差異表達基因數(shù)達16 407個(上調(diào)8 538，下調(diào)7 869)；QZT-QZCK間差異表達基因有16 021個(上調(diào)8 471，下調(diào)7 550)。MDT-QZT和MDCK-QZCK差異基因分別為6 027個(上調(diào)3 440，下調(diào)2 587)和5 994個(上調(diào)3 456，下調(diào)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個組間同步上調(diào)表達的差異基因為4 121個，同步下調(diào)表達的差異基因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個亞類，且均在分子功能過程中“細胞”、“膜”、細胞組分部分中“結(jié)合”和“催化活性”、生物學過程中“細胞過程”、“代謝過程”等功能中注釋基因相對較多，能被注釋功能的差異基因分別為13 524個和13 249個，MD較QZ多出2.08%。MDCK-MDT組在分子功能中“結(jié)合”、“催化活性”中差異基因有6 659個(上調(diào)3 669，下調(diào)2 990)和6 651個(上調(diào)3 306，下調(diào)3 345)。細胞組分中“細胞”、“膜”中差異基因分別為6 294個(上調(diào)3 357，下調(diào)2 937)和4 982個(上調(diào)2 281，下調(diào)2 701)。生物學過程中“細胞過程”和“代謝過程”注釋的基因相對較多，分別有5 079個(上調(diào)2 753，下調(diào)2 326)和4 334個(上調(diào)2 323，下調(diào)2 011)。前20條GO注釋如圖4 A所示。

QZCK-QZT組在生物學過程中，分子功能中“結(jié)合”、“催化活性”中有6 539個(上調(diào)3 608，下調(diào)2 931)和6 494個(上調(diào)3 234，下調(diào)3 260)基因。細胞組分中“細胞”、“膜”注釋的差異基因分別為6 176個(上調(diào)3 358，下調(diào)2 818)和4 877個(上調(diào)2 260，下調(diào)2 617)。生物學過程“細胞過程”、“代謝過程”富集較多，分別為5 001個(上調(diào)2 268，下調(diào)2 733)和4 309個(上調(diào)2 338，下調(diào)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. 結(jié)合；2. 催化活性；3. 轉(zhuǎn)運活性；4. 轉(zhuǎn)錄調(diào)節(jié)活性；5.結(jié)構(gòu)分子活性；6. 細胞；7膜;8細胞器；9. 膜部件；10. 細胞器部分；11. 含蛋白質(zhì)復合物；12. 膜腔；13. 細胞過程；14. 代謝過程；15.生物調(diào)節(jié)；16. 應(yīng)激響應(yīng)；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上調(diào)基因顯著富集GO分類條目194條，前2條為在細胞組成中“細胞”(GO:0005622)和生物學過程中“細胞過程”(GO:0043604)，下調(diào)基因顯著富集GO分類條目132條，前2條為細胞組成中“膜”(GO:0016020)和生物學過程中“細胞過程”(GO:0006464)。低溫脅迫誘導QZT-QZCK上調(diào)基因最顯著富集在生物過程中“細胞過程”(GO：0043604；GO:0006412)，下調(diào)基因顯著富集在細胞組成中“膜”(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個(上調(diào)1 573，下調(diào)1 435)，遺傳信息處理中“翻譯”1 064個(上調(diào)640，下調(diào)424)、有機系統(tǒng)中的“環(huán)境適應(yīng)”934個(上調(diào)503，下調(diào)431)、環(huán)境信息處理中的“信號轉(zhuǎn)導”643個(上調(diào)254，下調(diào)389)、細胞過程中“運輸和分解”592個(上調(diào)254，下調(diào)338)。如圖5A所示。

QZT-QZCK的DEGs分別為代謝過程中的“全球和概覽圖”2 978個(上調(diào)1 586，下調(diào)1 392)、遺傳信息處理中“翻譯”1 086個(上調(diào)665，下調(diào)421)、有機系統(tǒng)的“環(huán)境適應(yīng)”939個(上調(diào)486，下調(diào)453)，環(huán)境信息處理中的“信號轉(zhuǎn)導”651個(上調(diào)402，下調(diào)249)、細胞過程中“運輸和分解”573個(上調(diào)320，下調(diào)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)境適應(yīng);2.全球和概覽圖;3.碳水化合物代謝;4.脂質(zhì)代謝;5.其他次生代謝產(chǎn)物的生物合成;6.氨基酸代謝;7.能量代謝;8.輔助因子和維生素的代謝;9.其他氨基酸代謝;10.聚糖生物合成與代謝;11.萜類化合物和聚酮類化合物的代謝;12.核苷酸代謝;13.翻譯;14.折疊、分類和降解;15.轉(zhuǎn)錄;16.復制和修復;17.信號轉(zhuǎn)導;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個代謝通路，即核糖體、植物激素信號轉(zhuǎn)導通路、植物晝夜節(jié)律通路、光合作用天線蛋白通路、卟啉與葉綠素代謝通路，其中植物激素信號傳導通路最多，差異基因富集數(shù)量達543個，其次為核糖體通路為341個基因(表4)。QZT-QZCK組差異基因顯著富集在3個代謝途徑，即核糖體通路、磷酸肌醇代謝和磷脂酰肌醇信號系統(tǒng)，通路富集基因數(shù)分別為353，119和113個；富集差異基因最多的為植物激素信號傳導通路，達513個，其次為植物促分裂素原活化蛋白激酶(MAPK)信號通路，差異基因達387個，差異不顯著(表5)。耐低溫材料MD中顯著富集的KEGG通路較低溫敏感材料QZ中多，兩者均共同富集的通路為核糖體通路。植物激素信號轉(zhuǎn)導通路、植物晝夜節(jié)律通路、光合作用天線蛋白通路、卟啉與葉綠素代謝通路只在MD材料中顯著富集，說明這些通路可能參與耐寒材料響應(yīng)低溫脅迫的表達。

表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 耐低溫脅迫響應(yīng)的差異基因分析

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

2.7 其他與低溫脅迫相關(guān)的差異基因及通路

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

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

2.6 實時熒光定量PCR驗證

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

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

3 討論

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

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

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

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

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

4 結(jié)論

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