植原體侵染對桑樹葉片轉錄組的影響

楊金宏，孔衛青

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

0 引言

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

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

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

1 材料與方法

1.1 實驗材料

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

1.2 葉片RNA提取、cDNA文庫構建和轉錄組序列測定

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

1.3 測序數據處理與功能注釋

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

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

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

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

2 結果與分析

2.1 桑樹葉片轉錄組測序數據結果

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

表1 測序數據統計

2.2 差異表達基因分析

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

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

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

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

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

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

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

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

圖3 桑樹轉錄組DEGs的KEGG代謝途徑分類

2.5 抗病相關DEGs基因分析

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

表2 抗病相關差異表達基因數量統計

續表2 抗病相關差異表達基因數量統計

3 討論

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

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

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

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