楊金宏，孔衛(wèi)青

(安康學院陜西省蠶桑重點實驗室，陜西 安康 725000)

0 引言

桑樹黃化型萎縮病(Mulberry Yellow Dwarf Disease，MYDD)是桑樹的重大病害，染病桑樹枝短葉小，葉片異常黃化，病株逐漸衰竭而枯死，全國范圍內(nèi)的桑園均可受到該病的侵害，常導(dǎo)致大片桑園毀滅，嚴重制約蠶業(yè)生產(chǎn)的健康發(fā)展[1]。該病由植原體(Phytoplasma)引起，這是一類重要的原核致病菌，寄主范圍廣，能夠侵染98科1 000多種植物[2]，寄生于植物韌皮部，隨植物體內(nèi)營養(yǎng)流進行運動，先下行到根部，再在植物韌皮部篩管內(nèi)被動上行誘發(fā)癥狀出現(xiàn)，使病株葉片、葉脈和韌皮部部分細胞的細胞核及細胞器結(jié)構(gòu)發(fā)生改變，甚至解體[3]，危害嚴重[4,5]。

目前桑樹黃花型萎縮病的研究，主要集中在病株同工酶、代謝物、蛋白質(zhì)組和miRNA的變化的研究。雷國新等研究發(fā)現(xiàn)MYDD桑樹枝皮的過氧化物同工酶譜帶較正常植株有不同程度的加強，且隨染病后經(jīng)過日數(shù)的增多而愈發(fā)顯著[6]。這與過氧化物酶可參與植物的多種代謝途徑，被認為與植物抗病抗逆性關(guān)系密切一致。韓雪娟等利用質(zhì)譜(Mass Spectrometry，MS)分析了健康和MYDD植株葉片汁液中的極性代謝物，發(fā)現(xiàn)葉片中有65種代謝物的濃度發(fā)生了變化，大部分濃度降低[7]。袁傳忠分析桑樹韌皮部汁液中有56種代謝物的濃度發(fā)生變化[8]。冀憲領(lǐng)研究了桑樹葉片蛋白質(zhì)組響應(yīng)植原體侵染的變化，發(fā)現(xiàn)植原體侵染可以引起葉片中包括參與自由基清除、代謝、防御、碳水化合物和能量代謝等方面15種蛋白質(zhì)的變化[9]。蓋英萍等的研究顯示，健康和MYDD桑樹韌皮部汁液中miRNA的表達存在差異，且鑒定出部分桑樹韌皮部汁液響應(yīng)植原體侵染的miRNAs[10]。

植原體的侵染對桑樹代謝物、蛋白質(zhì)組、miRNAs表達等均有重要影響，植原體侵染與桑樹黃花型萎縮病的發(fā)生可能涉及到一系列基因的差異表達，影響到多種代謝過程，形成復(fù)雜的基因表達和代謝調(diào)控網(wǎng)絡(luò)。轉(zhuǎn)錄組(transcriptome)狹義上指某一生理條件下，細胞內(nèi)所有mRNA的集合，是連接遺傳信息與生物功能蛋白質(zhì)的必然紐帶，也是生物體最重要的基因調(diào)控方式?；诖?，本文對黃花型萎縮病桑樹葉片轉(zhuǎn)錄組的變化，及相應(yīng)的代謝途徑和基因進行研究，以揭示桑樹響應(yīng)植原體侵染的分子機理研究提供參考。

1 材料與方法

1.1 實驗材料

桑樹品種為湖桑32，選擇經(jīng)qPCR驗證的健康和有明顯感病性狀的桑樹各5株，用花盆種植于組培室中，分別取每株中部葉片3枚，每枚葉片切取0.1 g用液氮速凍，后保存于-80℃冰箱中備用。

1.2 葉片RNA提取、cDNA文庫構(gòu)建和轉(zhuǎn)錄組序列測定

均勻混合上述桑樹的葉片組織，液氮研磨，參照Invitrogen公司的TRIzol Reagent說明書提取桑樹葉片總RNA。使用1.0%瓊脂糖凝膠電泳和Nanodrop檢測RNA的純度(D260/D280)和濃度。經(jīng)檢測合格的RNA樣品進行文庫構(gòu)建。文庫構(gòu)建方法參照TruseqTMRNA sample prep Kit(Illumina公司)說明書進行，文庫構(gòu)建完成后，由上海美吉生物醫(yī)藥科技有限公司在IlluminaHiseq 2 000測序平臺完成2×100 bp的DNA雙末端測序。

1.3 測序數(shù)據(jù)處理與功能注釋

堿基讀取得到原始序列數(shù)據(jù)(raw data)后，對其進行質(zhì)控處理以獲得高質(zhì)量的測序結(jié)果(Clean Data)，包括去掉含有接頭Reads和低質(zhì)量Reads(質(zhì)量值Q小于20、含N較多和長度小于20nt)。后用Trinity軟件對Clean Data進行序列拼接和ORF預(yù)測，得到轉(zhuǎn)錄本和ORF序列特征和長度分布，BLAST軟件對預(yù)測的氨基酸序列和剩余未預(yù)測出ORF的序列與NCBI的非冗余蛋白數(shù)據(jù)庫、String、Gene數(shù)據(jù)庫進行比對(E-value：1e-5)，對預(yù)測ORF進行功能注釋。Bowtie軟件以Clean Data和川桑MorusnotabilisSchneid基因組為參考，獲得序列在參考基因組或基因上的位置信息及測序樣品特有的序列特征信息[11]。

后利用Blast2go軟件進行GO功能注釋，利用BLAST軟件與STRING 9.0數(shù)據(jù)庫(http://string-db.org)和KEGG genes數(shù)據(jù)庫(http://www.genome.jp/kegg/genes.html)比對，以所得COG(Clusters of Orthologous Groups of proteins)number進行功能歸類[12]，以得到的KO編號查找具體的生物學通路，獲得基因可能參與的所有生物學通路。最后使用RSEM軟件對UniGenes和基因的表達水平進行定量，以FPKM值作為基因表達水平指標。

1.4 差異表達基因篩選與顯著性富集分析

使用edgeR軟件進行差異表達基因篩選，符合錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)小于0.05且差異倍數(shù)(fold change)在3倍以上的基因被定義為樣品間差異表達基因(differential gene，DEG)[13]。進一步使用goatools軟件對對差異基因進行GO功能顯著性富集分析[14]，使用KOBAS軟件對差異表達基因進行通路富集分析，以尋找與研究背景密切關(guān)聯(lián)的生物學過程[15]，與植物抗性基因數(shù)據(jù)庫PRGDB(http://prgdb.org)比對[16]。

2 結(jié)果與分析

2.1 桑樹葉片轉(zhuǎn)錄組測序數(shù)據(jù)結(jié)果

對健康(MO)和病樹(MV)葉片轉(zhuǎn)錄組測序結(jié)果表明(表1)，2樣本分別獲得63 484 370和61 812 660條raw reads，質(zhì)控處理后，分別獲得2樣本的測序數(shù)據(jù)5.87Gb和5.70Gb數(shù)據(jù)，含61 230 412和59 551 764個reads片段，Q30堿基百分比分別為91.69%和91.55%，能夠滿足后續(xù)分析要求。

表1 測序數(shù)據(jù)統(tǒng)計

2.2 差異表達基因分析

根據(jù)FDR<0.05且Fold change>3的標準，以健康桑樹葉片為參照，篩選統(tǒng)計兩樣本間的差異表達基因(different expression genes，DEGs)，結(jié)果1 394個基因具有差異表達，其中1 032個表達上調(diào)，362個表達下調(diào)(圖1)，575個基因在MO樣本中有表達，1 250個在MV樣本中具有表達，兩樣本中均有表達基因數(shù)431個(圖1)。

圖1 兩樣本基因表達差異統(tǒng)計(FDR<0.05)

2.3 差異表達基因的GO功能分類

進一步分析兩樣品DEGs的GO功能分類分析，結(jié)果顯示(圖2)，分子功能(molecular function)分類中，主要聚集在結(jié)合(binding)和催化活性(catalytic activity)；細胞組分(cell component)分類主要聚集在細胞(cell)和膜組分(membrane part)；生物學過程(biological process)分類中主要集中在細胞進程(cellular process)和代謝過程(metabolic process)。

圖2 差異表達基因的GO注釋

2.4 差異表達基因的KEGG代謝通路分析

以KEGG數(shù)據(jù)庫作參考進行分析，結(jié)果其中的399個DEGs參與了代謝通路的6大類39個亞類175小類。圖3A為KEGG大類分區(qū)6大類分別占比例，可以看到，6大類中，代謝類通路所占比例最大，為43.30%，其次是人類疾病和遺傳信息處理相關(guān)通路，占比分別為16.29%和13.93%，細胞過程、生物系統(tǒng)和環(huán)境信息處理相關(guān)的通路，占比為8.25%、8.79%和9.43%。

圖3B是小類分區(qū)中存在大于10個DEGs的代謝通路的情況，其中小類中數(shù)量最多的是KO03010，來自遺傳信息處理大類的翻譯亞類的核糖體，總數(shù)達76個。KO00010KO00380來自代謝類通路，KO00010、KO00620、KO00640和KO00650為碳水化合物代謝亞類，KO00190、KO00680和KO00910來自能量代謝亞類，KO00280和KO00380來自氨基酸代謝亞類，這3個亞類分別占代謝通路的10.61%、7.07%和8.36%，是代謝類通路的主要組成部分。另外來自環(huán)境信息處理大類信號轉(zhuǎn)導(dǎo)亞類的KO02020數(shù)量也較多，為32個。

圖3 桑樹轉(zhuǎn)錄組DEGs的KEGG代謝途徑分類

2.5 抗病相關(guān)DEGs基因分析

植物抗病基因(R基因)起源于植物和動物共有的一種古老的免疫系統(tǒng)，在識別病原體特異性蛋白質(zhì)方面起著關(guān)鍵作用。進一步依據(jù)基因的功能注釋和已有的基因功能的報道，對可能與植物病害相關(guān)基因進行篩選，結(jié)果獲得89個基因(表2,續(xù)表2)，這些基因主要涉及絲氨酸/蘇氨酸蛋白激酶(Serine/theronine-protein kinase)、抗病蛋白(Putative disease resistance protein)、受體類蛋白激酶(receptor-like protein kinase)、脂質(zhì)氧化酶等37類。以正常未患病桑葉樣品為對照，結(jié)果顯示：抗病相關(guān)DEGs中，50個基因上調(diào)表達，39個基因下調(diào)表達。

表2 抗病相關(guān)差異表達基因數(shù)量統(tǒng)計

續(xù)表2 抗病相關(guān)差異表達基因數(shù)量統(tǒng)計

3 討論

桑樹基因組測序的完成為系統(tǒng)、深入地挖掘桑樹功能基因，揭示對病原應(yīng)答的分子機制提供了重要基礎(chǔ)條件[17]。在有參考基因組的基礎(chǔ)上進行轉(zhuǎn)錄組測序和數(shù)據(jù)分析，可以獲得更精準的序列數(shù)據(jù)信息。川?；蚪M序列共注釋基因29 338個[18]，筆者研究測序獲得了健康和MYDD桑樹葉片兩個樣本的11.57Gb數(shù)據(jù)，共得到62 953個UniGenes，注釋了27 464個基因，覆蓋了川?；蚪M注釋的83%的基因，進一步去除植物的組織特異表達，說明測序數(shù)據(jù)和組裝效果好，注釋信息豐富。將27 464個基因與COG數(shù)據(jù)庫識別并對直系同源基因產(chǎn)物進行分類，結(jié)果有21 844條與數(shù)據(jù)庫中的基因具有相似性，分別被注釋到23個COG功能分類，涉及了桑樹大部分的生命活動，與10～20%譜系特有基因結(jié)果一致[19]。同時發(fā)現(xiàn)了378個在川桑中沒有注釋的基因，而通過巴旦木、葡萄、可可、野草莓等植物的同源序列進行了注釋。這說明桑屬植物內(nèi)具有較高的遺傳多樣性，也為桑屬植物多樣性的研究和優(yōu)良育種基因的鑒定提供了遺傳背景參考。

為明確植原體侵染對桑樹葉片轉(zhuǎn)錄組的影響，本研究篩選到了一批在健康和MYDD桑樹葉片中存在差異表達的UniGenes，找到了2 766條和3 570條分別在健康和MYDD桑樹葉片中單獨表達的UniGenes。結(jié)合UniGenes的GO功能差異富集，找到了富集量超過50條UniGenes的41個GO分類條目，7 909條對應(yīng)關(guān)系，共涉及到11個亞類，發(fā)現(xiàn)植原體侵染相關(guān)的差異UniGenes參與了分子功能的結(jié)合類在其中占比最多，其次是生物學過程的細胞進程類和細胞組分的細胞等生物學過程，這與植原體侵染導(dǎo)致植物內(nèi)源激素失衡的研究有一定的對應(yīng)關(guān)系。植原體能夠造成寄主的內(nèi)源激素紊亂，細胞分裂素含量增加，赤霉素含量降低，玉米素/生長素比值增加，而內(nèi)源激素在調(diào)控植物細胞的分裂、伸長、分化、生長等生物學過程中起著重要作用[20]。

進一步對健康和MYDD病植株的差異表達UniGenes基因與生命和環(huán)境聯(lián)系起來，進行基于KEGG的基因和基因組的代謝通路富集分析，找到了8個差異顯著的通路，涉及40余個蛋白。由于基因是相互影響的，1個基因往往有多個功能，而具體功能往往是多個基因共同作用的結(jié)果[21]，所以富集分析的結(jié)果總是會涉及特別多的通路。而且植原體的侵染可以對寄主植物產(chǎn)生多方面的影響，如蛋白質(zhì)和氨基酸組成、IAA信號轉(zhuǎn)導(dǎo)、過氧化物酶、酯酶等酶系統(tǒng)，次生代謝的生物合成[22]。本研究的KO00910和KO00250途徑的GAD、GdhA等谷氨酸代謝相關(guān)基因與在甜櫻桃、棗中得到研究結(jié)果一致[23～24]，KO00908玉米素生物合成Zeatinbiosynthesi也在多個物種中得到一致的結(jié)果[24～25]。

Cardiac muscle contraction(KO04260)是一個在哺乳動物中注釋的細胞通路，這部分基因在植物中多是葉綠體和線粒體基因。葉綠體是植物光合作用的重要器官，在植物的生長發(fā)育中具有重要作用，許多植物在被植原體感染后，光合作用顯著降低，碳水化合物的積累被抑制[26]。本研究發(fā)現(xiàn)部分參與光合作用的基因在MYDD桑樹葉片中的表達出現(xiàn)下調(diào)，這可能造成電子傳遞受阻，引起感染植物的光合作用的降低[27]。