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

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

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

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

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

1 材料與方法

1.1 供試材料

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

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

1.2 轉錄組測序

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

1.3 功能注釋與分類

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

1.4 CDC相關基因的篩選和進化分析

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

2 結果與分析

2.1 “魯保一號”菌株轉錄組數(shù)據(jù)的組裝

基于邊合成邊測序（sequencing by synthesis，SBS）技術，使用Illumina Hiseq 2500高通量測序平臺對“魯保一號”菌株完成轉錄組測序工作。共獲得總長度為431 911 195 bp的序列信息，進一步組裝獲得10 013 398個contig序列，主要以長度為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）。

對所獲得的contig數(shù)據(jù)進行進一步組裝，得到總長度為61 674 287 bp的transcripts，共25 588條，其N50為4 038 bp，組裝完整性較高。長度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）。

對獲得的transcripts序列進行進一步組裝，得到17 031條unigenes序列，總長度為31 126 662 bp，平均長度為1 827.65 bp，N50長度為3 093 bp。長度為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序列與多個公共數(shù)據(jù)庫進行比對，其中KOG數(shù)據(jù)庫中共有5 228個unigene獲得注釋，GO數(shù)據(jù)庫中共5 192個，NR數(shù)據(jù)庫中共9 991個（表2）。共獲得10 538個有注釋信息的unigene序列，占全部unigene序列的61.9%。

表1 “魯保一號”菌株轉錄組數(shù)據(jù)的組裝統(tǒng)計

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

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

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

圖1 unigene的GO分類結果

圖2 “魯保一號”菌株的同源物種分布

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

圖3 unigene的KOG功能分類

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

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

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

3 討論與結論

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

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

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