基于高通量測(cè)序的菟絲子生防菌“魯保一號(hào)”轉(zhuǎn)錄組學(xué)研究

李健，李美，高興祥，房鋒

（山東省農(nóng)業(yè)科學(xué)院植物保護(hù)研究所，山東 濟(jì)南 250100）

菟絲子是菟絲子屬（Cuscuta sp.）植物的統(tǒng)稱，能夠通過吸器汲取寄主的養(yǎng)分，是危害嚴(yán)重的一類惡性寄生性雜草。近些年，菟絲子的侵入不僅對(duì)大豆等作物造成了危害，而且對(duì)我國(guó)甘肅、寧夏、內(nèi)蒙古和遼寧等地的草原牧區(qū)生態(tài)造成了嚴(yán)重威脅［1-3］。尖孢炭疽菌“魯保一號(hào)”［3］是山東省農(nóng)業(yè)科學(xué)院劉志海等人于1963年研制成功的、一種對(duì)菟絲子有特殊防效的微生物菌［4-6］。該菌曾在我國(guó)多個(gè)省、區(qū)得到推廣和應(yīng)用，對(duì)控制菟絲子的蔓延和危害起到了一定作用［4，5］。由于該菌株致病力易發(fā)生退化，且應(yīng)用較早，基因資源信息匱乏，嚴(yán)重影響了對(duì)其進(jìn)一步的研究和應(yīng)用［7］。因此，系統(tǒng)研究“魯保一號(hào)”菌株的轉(zhuǎn)錄組信息，挖掘相關(guān)功能基因，對(duì)解決其培養(yǎng)過程中的致病力退化問題至關(guān)重要。

轉(zhuǎn)錄組是指細(xì)胞在特定狀態(tài)下表達(dá)的全部RNA的總和，反映了相應(yīng)物種在一定狀態(tài)下的基因表達(dá)狀況。由轉(zhuǎn)錄組延伸而來的差異轉(zhuǎn)錄組能夠反映相同基因在不同條件下表達(dá)水平的差異，為揭示不同基因的相互調(diào)控模式及各自功能提供了可能［8，9］。隨著測(cè)序儀器的改進(jìn)和測(cè)序原理的明晰，目前新一代轉(zhuǎn)錄組測(cè)序（高通量測(cè)序，highthroughput sequencing）技術(shù)已經(jīng)得到了越來越廣泛的應(yīng)用［10，11］。該技術(shù)可以在芯片上并行對(duì)數(shù)百萬計(jì)的DNA分子進(jìn)行大規(guī)模測(cè)序，從而獲得海量的測(cè)序結(jié)果，結(jié)合相應(yīng)分析手段，使得對(duì)缺乏基因組信息的物種進(jìn)行細(xì)致全貌的分析成為可能，是對(duì)傳統(tǒng)測(cè)序技術(shù)的一次變革［12-14］。

“魯保一號(hào)”菌株具有較高的應(yīng)用價(jià)值，但是由于其開發(fā)時(shí)間較早，且受制于其菌種致病力退化，對(duì)該菌株的研究基本停止，而對(duì)其功能基因信息相關(guān)研究則未見報(bào)道。對(duì)于尚未系統(tǒng)開展基因組學(xué)研究的物種來說，獲得大通量的基因資源信息是解決問題的首要步驟，新一代高通量測(cè)序技術(shù)的發(fā)展為從分子生物學(xué)水平研究“魯保一號(hào)”菌株提供了便利。本試驗(yàn)首次將高通量測(cè)序技術(shù)應(yīng)用到“魯保一號(hào)”菌株轉(zhuǎn)錄組研究中，對(duì)測(cè)序獲得的海量數(shù)據(jù)進(jìn)行拼接與組裝，結(jié)合生物信息學(xué)技術(shù)對(duì)所獲得的unigene進(jìn)行功能注釋和功能分類；以測(cè)序獲得的數(shù)據(jù)信息為依據(jù)，篩選CDC（細(xì)胞分裂周期基因，cell division cycle gene）系列相關(guān)基因，并系統(tǒng)分析其遺傳進(jìn)化關(guān)系，為“魯保一號(hào)”菌株基因功能的研究奠定基礎(chǔ)。

1 材料與方法

1.1 供試材料

試驗(yàn)于2015年12月—2017年5月在山東省農(nóng)業(yè)科學(xué)院植物保護(hù)研究所雜草科學(xué)實(shí)驗(yàn)室開展。供試菌株為“魯保一號(hào)”，保存于山東省農(nóng)業(yè)科學(xué)院植物保護(hù)研究所雜草科學(xué)研究室。

培養(yǎng)基為馬鈴薯葡萄糖培養(yǎng)基（PDA）和完全培養(yǎng)基（CM）。

1.2 轉(zhuǎn)錄組測(cè)序

PDA培養(yǎng)基活化復(fù)壯后的“魯保一號(hào)”菌株記為L(zhǎng)B-1，連續(xù)繼代（培養(yǎng)皿內(nèi)生長(zhǎng)7 d為一代）培養(yǎng)10代后的致病力減弱菌株記為L(zhǎng)B-1-10。挑取菌落邊緣長(zhǎng)勢(shì)一致的菌絲，接種于液體CM培養(yǎng)基內(nèi)，160 r/min、黑暗搖培72 h，雙層紗布過濾收集菌絲，菌絲1∶1混合后經(jīng)液氮速凍，送至北京百邁克生物科技有限公司進(jìn)行RNA的提取與轉(zhuǎn)錄組學(xué)分析。獲得高質(zhì)量的原始測(cè)序數(shù)據(jù)后，通過Trinity組裝軟件［15］對(duì)相應(yīng)序列進(jìn)行組裝拼接。首先將測(cè)序reads打斷為較短的片段（kmer），然后將這些小片段通過序列拼接組裝成較長(zhǎng)的重疊群（contig），并利用這些片段之間的重疊，得到片段集合（component），最后利用de-Bruijn圖的方法和測(cè)序read信息，在各個(gè)片段集合中分別識(shí)別轉(zhuǎn)錄本（transcript）序列，對(duì)轉(zhuǎn)錄本進(jìn)行同源聚類和拼接得到單基因簇（unigene）。

1.3 功能注釋與分類

使用BLAST軟件［16］將測(cè)序獲得的unigene序列與NCBI的非冗余核酸數(shù)據(jù)庫（non-redundant protein database，NR）、Swiss-Prot（swissprot protein sequence database）和蛋白質(zhì)直系同源數(shù)據(jù)庫（cluster of orthologous groups，COG）［17-19］等蛋白質(zhì)數(shù)據(jù)庫進(jìn)行比對(duì)分析，獲得最佳功能注釋。之后使用HMMER軟件［20］與Pfam［21］數(shù)據(jù)庫比對(duì)，獲得unigene的注釋信息。根據(jù)NCBI數(shù)據(jù)庫的功能注釋信息，使用GO軟件［22］得到unigene的GO條目，然后用WEGO軟件［23］進(jìn)行分類統(tǒng)計(jì)。數(shù)據(jù)分析過程中，選擇BLAST參數(shù)E-value不大于1e-5和HMMER參數(shù)E-value不大于1e-10。

1.4 CDC相關(guān)基因的篩選和進(jìn)化分析

前期研究顯示，細(xì)胞分裂紊亂是“魯保一號(hào)”菌株連續(xù)繼代培養(yǎng)后的重要現(xiàn)象［7］。根據(jù)對(duì)轉(zhuǎn)錄組數(shù)據(jù)的功能注釋結(jié)果，結(jié)合比對(duì)分析篩選獲得了6個(gè)CDC相關(guān)基因。利用MEGA 6軟件對(duì)獲得的相關(guān)基因進(jìn)行系統(tǒng)進(jìn)化樹分析。

2 結(jié)果與分析

2.1 “魯保一號(hào)”菌株轉(zhuǎn)錄組數(shù)據(jù)的組裝

基于邊合成邊測(cè)序（sequencing by synthesis，SBS）技術(shù)，使用Illumina Hiseq 2500高通量測(cè)序平臺(tái)對(duì)“魯保一號(hào)”菌株完成轉(zhuǎn)錄組測(cè)序工作。共獲得總長(zhǎng)度為431 911 195 bp的序列信息，進(jìn)一步組裝獲得10 013 398個(gè)contig序列，主要以長(zhǎng)度為200～300 bp的contig序列為主，有9 962 103條，占總體的99.49%；300～500 bp的contig序列有22 853條，占總體的0.23%；500～1 000 bp的contig序列有14 422條，占總體的0.14%；1 000～2 000 bp的contig序列有7 677條，占總體的0.08%；≥2 000 bp的contig序列有6 343條，占總體的0.06%（表1）。

對(duì)所獲得的contig數(shù)據(jù)進(jìn)行進(jìn)一步組裝，得到總長(zhǎng)度為61 674 287 bp的transcripts，共25 588條，其N50為4 038 bp，組裝完整性較高。長(zhǎng)度200～300、300～500、500～1 000、1 000～2 000 bp和≥2 000 bp的transcripts序列分別占總體的9.27%、10.70%、15.82%、21.90%和42.30%（表1）。

對(duì)獲得的transcripts序列進(jìn)行進(jìn)一步組裝，得到17 031條unigenes序列，總長(zhǎng)度為31 126 662 bp，平均長(zhǎng)度為1 827.65 bp，N50長(zhǎng)度為3 093 bp。長(zhǎng)度為200～300、300～500、500～1 000、1 000～2 000 bp和≥2 000 bp的unigene序列分別占總體的12.87%、13.47%、18.22%、22.99%和32.45%（表1）。

2.2 unigene的功能注釋、分類和代謝途徑分析

將拼裝得到的unigene序列與多個(gè)公共數(shù)據(jù)庫進(jìn)行比對(duì)，其中KOG數(shù)據(jù)庫中共有5 228個(gè)unigene獲得注釋，GO數(shù)據(jù)庫中共5 192個(gè)，NR數(shù)據(jù)庫中共9 991個(gè)（表2）。共獲得10 538個(gè)有注釋信息的unigene序列，占全部unigene序列的61.9%。

表1 “魯保一號(hào)”菌株轉(zhuǎn)錄組數(shù)據(jù)的組裝統(tǒng)計(jì)

表2 BLAST比對(duì)公共數(shù)據(jù)庫結(jié)果

2.2.1 unigene的GO注釋結(jié)果 GO（基因本體，gene ontology）是一個(gè)被廣泛應(yīng)用的標(biāo)準(zhǔn)化基因功能分類數(shù)據(jù)庫，數(shù)據(jù)庫分類注釋結(jié)果總共有三大類，分別是分子功能（molecular function）、細(xì)胞組分（cellular component）和生物學(xué)過程（biological process）。本試驗(yàn)結(jié)果表明，可將轉(zhuǎn)錄組獲得的所有unigene劃分為52個(gè)功能組，其中3 754個(gè)屬于細(xì)胞組分，7 935個(gè)屬于分子功能，7 410個(gè)屬于生物學(xué)過程。其中細(xì)胞成分、細(xì)胞器成分、催化活性、結(jié)合活性、代謝進(jìn)程、細(xì)胞進(jìn)程和單一生物進(jìn)程涉及的unigene較多，而病毒體、胞外基質(zhì)、金屬伴侶蛋白活性、通道調(diào)節(jié)活性和細(xì)胞殺傷等涉及的unigene沒有或極少（圖1）。

2.2.2 unigene的NR注釋結(jié)果 使用BLAST軟件將unigene序列與NR數(shù)據(jù)庫比對(duì)，進(jìn)行序列相似性分析，得到與給定unigene具有最高序列相似性的蛋白描述，并獲得unigene蛋白的功能注釋信息。由圖2可知，86.04%的序列與已知炭疽菌序列有不同程度的同源性，相似序列匹配的近緣物種還有大豆、高粱、西瓜等，其他物種占13.05%。

圖1 unigene的GO分類結(jié)果

圖2 “魯保一號(hào)”菌株的同源物種分布

2.2.3 unigene的KOG注釋結(jié)果 “魯保一號(hào)”菌株的unigene根據(jù)其功能大致分為25類（圖3），涉及了大多數(shù)生命活動(dòng)，如RNA加工與修飾，染色體結(jié)構(gòu)和動(dòng)力學(xué)，能量產(chǎn)生與運(yùn)輸，細(xì)胞周期控制、細(xì)胞分裂及染色體分裂，氨基酸運(yùn)輸及代謝等。其中注釋最多的是一般功能預(yù)測(cè)類基因（R），其次是翻譯后修飾、蛋白折疊和分子伴侶類基因（O），再者是翻譯、核糖體結(jié)構(gòu)和生物發(fā)生類基因（J），只有極少數(shù)的細(xì)胞活性類基因（N）和胞外結(jié)構(gòu)類基因（W）。

圖3 unigene的KOG功能分類

2.3 CDC基因系統(tǒng)進(jìn)化樹分析

系統(tǒng)進(jìn)化樹分析表明（圖4），篩選獲得的6個(gè)CDC相關(guān)基因被分為三個(gè)亞組。其中CDC3和CDC6在同一亞組，CDC1、CDC5和CDC2在同一亞組，CDC4單獨(dú)一個(gè)亞組。

圖4 “魯保一號(hào)”菌株CDC基因系統(tǒng)進(jìn)化樹分析

3 討論與結(jié)論

菟絲子能夠通過吸器寄生多種作物，嚴(yán)重影響大豆（Glycine max）、牧草和蔬菜的產(chǎn)量和品質(zhì)［1-3］。“魯保一號(hào)”菌株對(duì)菟絲子防效良好，但在經(jīng)過一段時(shí)間的應(yīng)用后，受制于其致病力退化問題，最終被遺棄［6，7］。由于該菌發(fā)現(xiàn)、應(yīng)用較早，受制于當(dāng)時(shí)的技術(shù)條件，并未得到深入的遺傳學(xué)研究，也沒有關(guān)于其功能基因研究的報(bào)道［3，6，7］。開展該菌株的轉(zhuǎn)錄組學(xué)研究，初步獲得其轉(zhuǎn)錄組信息，對(duì)于挖據(jù)優(yōu)良基因資源、解決其致病力退化問題等具有重要意義。

轉(zhuǎn)錄組測(cè)序中獲得的unigene片段太短會(huì)導(dǎo)致在后續(xù)比對(duì)注釋過程中無法找到匹配序列。本研究采用Illumina Hiseq 2500高通量測(cè)序技術(shù)首次對(duì)“魯保一號(hào)”菌株的轉(zhuǎn)錄組進(jìn)行測(cè)序和組裝分析，共獲得17 031條unigene序列信息，平均長(zhǎng)度為1 827.65 bp，能夠很好地用于進(jìn)一步的功能注釋和分析，為后續(xù)批量分析“魯保一號(hào)”菌株功能基因提供了可能。結(jié)合相關(guān)生物信息學(xué)分析方法，對(duì)獲得的“魯保一號(hào)”菌株unigene序列信息與各數(shù)據(jù)庫進(jìn)行比對(duì)，進(jìn)行序列相似性和功能注釋分析。NR數(shù)據(jù)庫比對(duì)顯示，86.04%unigene標(biāo)注信息與炭疽菌序列一致，這也進(jìn)一步證明了“魯保一號(hào)”菌株為炭疽菌屬；另一方面，13.96%的unigene與其它物種有不同程度的同源性，為“魯保一號(hào)”菌株功能基因的進(jìn)一步挖掘提供了參考。GO分類進(jìn)一步顯示了“魯保一號(hào)”菌株生長(zhǎng)發(fā)育過程中基因表達(dá)譜的總體情況，其中細(xì)胞組分中的細(xì)胞成分、細(xì)胞器成分，分子功能中的催化活性、結(jié)合活性，生物學(xué)進(jìn)程中的代謝進(jìn)程、細(xì)胞進(jìn)程和單一生物進(jìn)程涉及的unigene較多，為下一步大量挖掘相關(guān)功能基因奠定了基礎(chǔ)?；赟SR位點(diǎn)的分子標(biāo)記在物種遺傳圖譜構(gòu)建、遺傳多樣性分析、相關(guān)生物進(jìn)程功能基因發(fā)現(xiàn)和以SSR為分子標(biāo)記的輔助育種等研究中得到了較為廣泛的應(yīng)用［24-26］。本研究通過查找發(fā)現(xiàn)了3 587個(gè)SSR位點(diǎn)，接下來可設(shè)計(jì)并篩選SSR引物，為進(jìn)一步開發(fā)新的SSR標(biāo)記奠定基礎(chǔ)。

前期報(bào)道顯示，“魯保一號(hào)”菌株連續(xù)培養(yǎng)后存在細(xì)胞分裂異常現(xiàn)象［7］。本研究以獲得的轉(zhuǎn)錄組數(shù)據(jù)庫為基礎(chǔ)，經(jīng)過序列比對(duì)分析，初步篩選獲得了6個(gè)CDC相關(guān)基因，分屬3個(gè)亞組。絲狀真菌內(nèi)CDC基因的功能研究表明，不同的CDC基因在生物功能發(fā)揮過程中起到不同作用［27］；同一亞組內(nèi)的基因在功能上存在相近或互補(bǔ)的可能，這為下一步集中開展“魯保一號(hào)”菌株CDC基因功能分析提供了參照。實(shí)時(shí)定量分析表明其中的兩個(gè)基因在連續(xù)培養(yǎng)后表達(dá)量顯著增加，說明這兩個(gè)基因可能與連續(xù)培養(yǎng)后的細(xì)胞分裂異常相關(guān)，可能參與了“魯保一號(hào)”菌株致病力退化的調(diào)控過程，為下一步的深入研究確立了目標(biāo)。本研究通過高通量測(cè)序獲得了“魯保一號(hào)”菌株的大量轉(zhuǎn)錄組信息，為其基因克隆、分子標(biāo)記發(fā)掘和基因組學(xué)研究等提供了有價(jià)值的數(shù)據(jù)。