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

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

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

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

1 材料與方法

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

供試小麥品種為晉麥47，其對低氮條件具有較好適應(yīng)性[17]，由山西省臨汾市農(nóng)業(yè)種子管理站提供。小麥種子經(jīng)75%乙醇表面消毒，并用蒸餾水沖洗4～5次，置于濕潤濾紙在24 ℃黑暗條件下培養(yǎng)至萌發(fā)，然后轉(zhuǎn)至人工氣候箱光照培養(yǎng)(16 h光照/8 h黑暗)，生長至二葉一心時選取整齊一致的幼苗，分別轉(zhuǎn)移至限氮營養(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濃度調(diào)整為0.05 mmol·L-1，其他成分與氮充分營養(yǎng)液相同。為保持小麥生長條件恒定，每3 d換一次營養(yǎng)液。

1.2 轉(zhuǎn)錄組測序

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

1.3 轉(zhuǎn)錄組數(shù)據(jù)的組裝

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

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

1.4 差異表達(dá)基因分析

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

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

在GO(gene ontology)數(shù)據(jù)庫對差異表達(dá)基因進(jìn)行功能注釋和富集分析。利用KEGG(kyoto encyclopedia of genes and genomes)對差異表達(dá)基因進(jìn)行通路注釋分析。應(yīng)用STRING軟件對上調(diào)和下調(diào)的差異表達(dá)基因進(jìn)行了蛋白質(zhì)相互作用分析。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)分析

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

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

2.2 轉(zhuǎn)錄組NR庫注釋

將獲得的轉(zhuǎn)錄組在NR數(shù)據(jù)庫中進(jìn)行比對，發(fā)現(xiàn)有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%。根據(jù)Blastx的結(jié)果可知，98%轉(zhuǎn)錄組序列的相似度大于60%，84%轉(zhuǎn)錄組序列的相似度大于80%。

2.3 轉(zhuǎn)錄組COG注釋

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

2.4 低氮脅迫下差異表達(dá)基因的篩選

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

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

2.5 葉片差異表達(dá)基因GO富集分析

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

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

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

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

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

注：A—RNA 加工與修飾；B—染色質(zhì)的結(jié)構(gòu)和動力學(xué)；C—能量產(chǎn)生和轉(zhuǎn)化；D—細(xì)胞周期控制、細(xì)胞分裂和染色體分區(qū)；E—氨基酸的運輸和代謝；F—核苷酸的運輸和代謝；G—碳水化合物的運輸和代謝；H—輔酶的運輸和代謝；I—脂質(zhì)運輸和代謝；J—翻譯，核糖體結(jié)構(gòu)和生物合成；K—轉(zhuǎn)錄；L—復(fù)制，重組和修復(fù)；M—細(xì)胞壁/膜/包膜的生物發(fā)生；N—細(xì)胞運動；O—翻譯后修飾，蛋白質(zhì)更新，分子伴侶；P—無機(jī)離子的運輸和代謝；Q—次生代謝產(chǎn)物的生物合成，轉(zhuǎn)運和分解代謝；R—一般功能預(yù)測；S—功能未知；T—信號轉(zhuǎn)導(dǎo)機(jī)制；U—細(xì)胞內(nèi)運輸，分泌和囊泡運輸；V—防御機(jī)制；W—細(xì)胞外結(jié)構(gòu)；Y—核結(jié)構(gòu)；Z—細(xì)胞骨架。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 氮脅迫條件下小麥葉片轉(zhuǎn)錄組COG注釋Fig.1 Function classification in clusters of orthologous groups (COG) of proteins

2.7 轉(zhuǎn)錄因子

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

表3 轉(zhuǎn)錄因子差異表達(dá)基因Table 3 Differentially expressed genes of transcription factor

2.8 差異表達(dá)基因蛋白互作網(wǎng)絡(luò)分析

蛋白質(zhì)相互作用在生命活動中起核心作用，由蛋白質(zhì)相互作用構(gòu)成的PPI網(wǎng)絡(luò)的拓?fù)涮匦苑治鍪呛蠡蚪M時代最重要內(nèi)容。對上調(diào)和下調(diào)的差異表達(dá)基因進(jìn)行了蛋白質(zhì)相互作用分析，在葉片的上調(diào)差異表達(dá)基因中未找到蛋白質(zhì)相互作用。如圖3所示，下調(diào)差異表達(dá)基因的PPI網(wǎng)絡(luò)分析中，BRADI2G40600.1與BRADI1G77050.1兩個結(jié)點分別與14個、12個蛋白相結(jié)合。進(jìn)一步分析表明BRADI2G40600.1結(jié)合點參與6個GO功能，分別為DNA修復(fù)(GO: 0006260)、核苷酸綁結(jié)合(GO: 0000166)、DNA結(jié)合(GO: 0003677)、ATP結(jié)合(GO: 0005524)、ATP酶活性(GO: 0016887)和核苷酸三磷酸酶活性(GO: 0017111)。

圖3 差異表達(dá)基因的蛋白質(zhì)互作網(wǎng)絡(luò)Fig.3 Protein interaction network of differentially expressed genes

3 討論

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

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