氮脅迫下小麥葉片轉錄組分析

張文云, 張建誠, 姚景珍*

(1.臨汾市農業(yè)種子管理站, 山西 臨汾 041000; 2.山西農業(yè)大學棉花研究所, 山西 運城 044000)

我國是世界小麥種植面積較大、總產量較高的國家之一，近年來種植面積已達2 427萬hm2，總產量約為1.2億t，分別約占世界的13%和19%[1]。氮肥作為提高作物產量和改善品質的主要栽培手段，在農業(yè)生產中起著重要作用[2]。在小麥生產中，為了追求高產往往過量施用氮肥，造成氮肥利用率降低，不僅導致資源浪費，而且給生態(tài)環(huán)境帶來巨大壓力。因此，深入研究作物體內的氮素循環(huán)和提高作物氮素利用效率的技術已成為當前國際的研究熱點。

轉錄組(transcriptome)是某個物種或者特定細胞類型產生的所有轉錄產物的集合，包括信使RNAs、非編碼RNAs和小RNA[11]。轉錄組學(tanscriptomics)是功能基因組學的重要組成部分，從整體水平上研究細胞中基因轉錄情況及轉錄調控規(guī)律，對特定條件下基因表達水平改變做定量分析[12]。轉錄組學獲得的大量數據經過專業(yè)的生物信息學分析，既可以還原細胞內基因表達的特征，也可以獲得多種基因表達調控的重要信息。目前，轉錄組學分析已應用于小麥在不同脅迫條件的響應[13-14]，孟晨[15]對耐鹽堿小麥山融4號進行轉錄組分析，發(fā)現在堿脅迫條件下其根中有大量被誘導表達的基因，部分堿脅迫誘導表達基因提高了轉基因植物的耐鹽堿性。Li等[16]利用轉錄組測序分析田間小麥響應冷脅迫的分子機制，發(fā)現CBF和ABA通路對小麥的冷適應非常重要，且在小麥6A和6D染色體上與冷脅迫有關的差異表達基因分布最多。目前通過轉錄組測序研究小麥適應低氮脅迫的分子機制的相關報道較少。因此，本文對氮脅迫下的小麥葉片轉錄組進行了研究和分析，旨在闡明小麥在氮脅迫下的基因轉錄和調控規(guī)律。

1 材料與方法

1.1 材料培養(yǎng)及處理

供試小麥品種為晉麥47，其對低氮條件具有較好適應性[17]，由山西省臨汾市農業(yè)種子管理站提供。小麥種子經75%乙醇表面消毒，并用蒸餾水沖洗4～5次，置于濕潤濾紙在24 ℃黑暗條件下培養(yǎng)至萌發(fā)，然后轉至人工氣候箱光照培養(yǎng)(16 h光照/8 h黑暗)，生長至二葉一心時選取整齊一致的幼苗，分別轉移至限氮營養(yǎng)液和氮充分營養(yǎng)液中生長。

氮充分營養(yǎng)液：0.2 mmol·L-1KH2PO4、1.0 mmol·L-1MgSO4·7H2O、0.75 mmol·L-1K2SO4、2.5 mmol·L-1CaCl2、0.000 05 mmol·L-1(NH4)6Mo7O24·4H2O、0.001 mmol·L-1H3BO3、0.000 5 mmol·L-1CuSO4·5H2O、0.001 mmol·L-1ZnSO4·7H2O、0.001 mmol·L-1MnSO4·H2O、0.1 mmol·L-1Fe-EDTA、1.0 mmol·L-1NH4NO3，pH 5.5。限氮營養(yǎng)液把NH4NO3濃度調整為0.05 mmol·L-1，其他成分與氮充分營養(yǎng)液相同。為保持小麥生長條件恒定，每3 d換一次營養(yǎng)液。

1.2 轉錄組測序

待小麥幼苗在限氮營養(yǎng)液中和氮充分營養(yǎng)液生長2周后取樣，應用天為時代公司RNA提取試劑盒提取總RNA。由北京百邁客生物科技有限公司采用Illumina-Hiseq測序平臺進行轉錄組測序。

1.3 轉錄組數據的組裝

為保證后續(xù)組裝質量，提高生物信息分析結果的科學性，應用SeqPrep和Sickle軟件去除原始測序數據中的測序接頭、引物序列和低質量值的讀段，處理后得到高質量測序數據。

然后，對原始數據過濾得到的RNA-seq測序數據進行從頭組裝。使用軟件為Trinity、Inchworm搜索方法建立K-mer graph (K=25)中全部可能路徑；采用Chrysalis對上述得到可能路徑進行分組后，對Chrysalis生成路徑進行修剪和整合，保存主要路徑。

1.4 差異表達基因分析

將測序得到的Reads與Unigene庫進行比對，根據比對結果結合RSEM進行表達量水平估計。利用RPKM (reads per kilobase per million reads)表示對應Unigene的表達豐度[18]。

在差異表達分析中采用校正后的P值，即FDR(false discovery rate)<0.05,且差異倍數FC(fold change)≥2作為差異表達基因(DEGs)的篩選標準。

在GO(gene ontology)數據庫對差異表達基因進行功能注釋和富集分析。利用KEGG(kyoto encyclopedia of genes and genomes)對差異表達基因進行通路注釋分析。應用STRING軟件對上調和下調的差異表達基因進行了蛋白質相互作用分析。

2 結果與分析

2.1 測序數據分析

利用Illumina HiseqTM2500高通量測序，從氮充分營養(yǎng)液生長的小麥葉片(CKL)和限氮營養(yǎng)液生長的小麥(NTL)樣本中分別獲得53 920 024、53 817 412條原始序列，平均讀長100 bp，原始堿基數分別為5 445 922 424、5 435 558 612 bp，堿基質量分值達到20分(測序錯誤率小于1%)的堿基比例均超過95%。

將低質量數據過濾后，從CKL和NTL樣本中獲得高質量序列分別為52 383 726和52 192 061條，高質量堿基數分別為5 108 134 511、5 091 803 373 bp，堿基質量分值達到20分(測序錯誤率小于1%)的堿基比例均超過99%。

2.2 轉錄組NR庫注釋

將獲得的轉錄組在NR數據庫中進行比對，發(fā)現有21%的isogene注釋到粗山羊草中，15.1%注釋到大麥中，13.6%注釋到烏拉爾圖小麥中，其余分別為二穗短柄草5.4%、水稻2.5%、小麥2.0%、高粱0.77%和玉米0.61%。15%的E值分布在10-20～10-5之間，6%的E值位于10-30～10-20之間，小于10-30的E值占78%。根據Blastx的結果可知，98%轉錄組序列的相似度大于60%，84%轉錄組序列的相似度大于80%。

2.3 轉錄組COG注釋

如圖1所示，基于序列同源性，27 194個isogene得到了COG注釋。所參與的代謝過程共25個分類，主要包括：一般功能預測，復制、重組和修飾，轉錄，信號傳遞機制，翻譯后修飾、蛋白翻轉、伴侶，翻譯、核糖體結構和生物發(fā)生，碳水化合物轉運和代謝，氨基酸轉運和代謝等。

2.4 低氮脅迫下差異表達基因的篩選

對轉錄組的表達量進行計算，并篩選差異表達基因(DEGs)后發(fā)現，低氮脅迫下葉片差異表達基因為1 267個，其中179個基因上調表達，1 088個基因下調表達(圖2)。

圖2 低氮脅迫下葉片差異表達基因散點圖Fig.2 Scatter plot of differentially expressed genes under low nitrogen stress

2.5 葉片差異表達基因GO富集分析

對氮脅迫條件下小麥差異表達基因進一步進行GO富集分析，確定差異表達基因主要行使的生物學功能。葉片中差異表達基因在GO注釋中共分為3個功能組(表1)：生物過程、細胞組分和分子功能。生物過程中，注釋分類到與植物光合作用相關過程的差異表達基因較多；細胞組分中，差異表達基因主要分布在質體、葉綠體類囊體膜(腔)等部位。以上結果表明，低氮脅迫會導致小麥葉片光合作用等過程發(fā)生顯著變化。

表1 低氮脅迫下葉片差異表達基因的GO功能注釋Table 1 GO functional annotation of differentially expressed genes under low nitrogen stress

2.6 葉片差異表達基因KEGG通路顯著性富集分析

葉片差異表達基因被注釋到 KEGG 數據庫6個一級通路的38個二級通路中，更進一步可細分到178個三級代謝或信號轉導途徑中(表2)。其中二級通路中氨基酸代謝(amino acid metabolism)、碳水化合物代謝(carbohydrate metabolism)、脂類代謝(lipid metabolism)、信號傳導(signal transduction)、免疫系統(tǒng)(immune system)被注釋的三級通路數目較多。氨基酸代謝、碳水化合物代謝和脂類代謝是植物生物合成的主要成員，其三級通路數目均為10左右，表明氮脅迫對這些通路影響較大，導致小麥的營養(yǎng)生長受到抑制。信號轉導的三級通路數目為15，則說明氮代謝在三級通路水平上開始進行級聯(lián)的響應，通過信號轉導影響更多的代謝過程。

表2 葉片差異表達基因KEGG通路顯著性富集分析Table 2 Enriched KEGG pathways of differentially expressed genes in leaf

注：A—RNA 加工與修飾；B—染色質的結構和動力學；C—能量產生和轉化；D—細胞周期控制、細胞分裂和染色體分區(qū)；E—氨基酸的運輸和代謝；F—核苷酸的運輸和代謝；G—碳水化合物的運輸和代謝；H—輔酶的運輸和代謝；I—脂質運輸和代謝；J—翻譯，核糖體結構和生物合成；K—轉錄；L—復制，重組和修復；M—細胞壁/膜/包膜的生物發(fā)生；N—細胞運動；O—翻譯后修飾，蛋白質更新，分子伴侶；P—無機離子的運輸和代謝；Q—次生代謝產物的生物合成，轉運和分解代謝；R—一般功能預測；S—功能未知；T—信號轉導機制；U—細胞內運輸，分泌和囊泡運輸；V—防御機制；W—細胞外結構；Y—核結構；Z—細胞骨架。Note: A—RNA processing and modification; B—Chromatin structure and dynamics; C—Energy production and conversion; D—Cell cycle control, cell division, chromosome partitioning; E—Amino acid transport and metabolism; F—Nucleotide transport and metabolism; G—Carbohydrate transport and metabolism; H—Coenzyme transport and metabolism; I—Lipid transport and metabolism; J—Translation, ribosomal structure and biogenesis; K—Transcription; L—Replication, recombination and repair; M—Cell wall/membrane/envelope biogenesis; N—Cell motility; O—Posttranslational modification, protein turnover, chaperones; P—Inorganic ion transport and metabolism; Q—Secondary metabolites biosynthesis, transport and catabolism; R—General function prediction only; S—Function unknown; T—Signal transduction mechanisms; U—Intracellular trafficking, secretion, and vesicular transport; V—Defense mechanisms; W—Extracellular structures; Y—Nuclear structure; Z—Cytoskeleton.圖1 氮脅迫條件下小麥葉片轉錄組COG注釋Fig.1 Function classification in clusters of orthologous groups (COG) of proteins

2.7 轉錄因子

大量研究表明，轉錄因子在植物生長或發(fā)育的各個階段以及適應各種逆境脅迫過程中都發(fā)揮著重要作用。低氮脅迫條件下，鑒定出轉錄因子類型有WRKY轉錄因子、MYB轉錄因子、乙烯反應轉錄因子和NAC等轉錄因子。經分析后發(fā)現，發(fā)生差異表達的共有23個，其中4個轉錄因子上調表達，19個下調表達(表3)。這些轉錄因子通過影響下游氮轉運相關基因的表達，參與植物氮缺乏脅迫的應答調控[19-21]。

表3 轉錄因子差異表達基因Table 3 Differentially expressed genes of transcription factor

2.8 差異表達基因蛋白互作網絡分析

蛋白質相互作用在生命活動中起核心作用，由蛋白質相互作用構成的PPI網絡的拓撲特性分析是后基因組時代最重要內容。對上調和下調的差異表達基因進行了蛋白質相互作用分析，在葉片的上調差異表達基因中未找到蛋白質相互作用。如圖3所示，下調差異表達基因的PPI網絡分析中，BRADI2G40600.1與BRADI1G77050.1兩個結點分別與14個、12個蛋白相結合。進一步分析表明BRADI2G40600.1結合點參與6個GO功能，分別為DNA修復(GO: 0006260)、核苷酸綁結合(GO: 0000166)、DNA結合(GO: 0003677)、ATP結合(GO: 0005524)、ATP酶活性(GO: 0016887)和核苷酸三磷酸酶活性(GO: 0017111)。

圖3 差異表達基因的蛋白質互作網絡Fig.3 Protein interaction network of differentially expressed genes

3 討論

比較不同生長環(huán)境、組織、器官和發(fā)育階段樣本間的基因表達差異是揭示特定分子機制的有效方法。本文對氮脅迫下轉錄組的表達量進行統(tǒng)計并標準化處理、篩選后，發(fā)現低氮脅迫下葉片差異表達基因為1 267個，其中葉片中179個基因上調表達，1 088個基因下調表達。練興明等[22]利用cDNA芯片技術分析了水稻幼苗中10 422個基因在受到低氮脅迫后的反應，3個時間點上共有471個(4.5%)基因表達量發(fā)生了改變，其中上調表達的有158個，下調表達的有355個。蔡紅梅等[23]在研究水稻低氮脅迫時發(fā)現，所檢測的32 341個基因中，有3 518個(10.88%)發(fā)生了差異表達，其中根中有2 214個，莖中有1 766個基因，根和莖共同基因有462個。本研究所檢測的27 194個基因中1 267個(4.66%)發(fā)生差異表達，差異表達基因比例與上述研究結果近似，低氮脅迫的共同反應相同，但數目低于后者研究結果，這些基因可能是導致晉麥47氮高效利用的主要原因。

PPI網絡分析表明BRADI2G40600.1 和 BRADI1G77050.1是氮脅迫下小麥葉片下調差異表達基因的結點。進一步分析表明BRADI1G77050.1與鈣離子結合相關。鈣離子是第二信使，參與真核生物對生物和非生物脅迫信號傳導途徑[28]。鈣參與植物脅迫反應中各種生理過程調節(jié)，如水和溶質的運動，細胞分裂，細胞壁的合成，呼吸和置換[29]。研究結果表明，BRADI1G77050.1通過鈣離子參與小麥氮脅迫反應。另外，BRADI2G40600.1在PPI網絡中與14個蛋白相關聯(lián)并且參與DNA復制，其在細胞循環(huán)和凋亡、植物生長發(fā)育和代謝過程中發(fā)揮重要作用[11,30]。