轉(zhuǎn)Cry1Ab/Cry1Ac基因大米粉飼喂印度谷螟后的比較轉(zhuǎn)錄組分析

唐培安，陶冶心，王康旭，吳學(xué)友，陳二虎

（南京財(cái)經(jīng)大學(xué) 食品科學(xué)與工程學(xué)院，江蘇高校糧油質(zhì)量安全控制及深加工重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023）

印度谷螟Plodia interpunctella(Hübener)是一種全球范圍內(nèi)廣泛發(fā)生的儲糧害蟲，不僅能夠通過直接取食對糧食造成損失，其生命活動過程還可以提高儲糧溫度，從而引起糧堆霉變和品質(zhì)劣變。蘇云金芽孢桿菌Bacillus thuringiensis（Bt）是一種革蘭氏陽性細(xì)菌，能夠產(chǎn)生具有殺蟲效果的晶狀蛋白（cry），對鱗翅目和鞘翅目害蟲具有明顯的防治效果[1]。隨著轉(zhuǎn)基因技術(shù)的不斷發(fā)展，利用轉(zhuǎn)基因進(jìn)行 cry蛋白異源表達(dá)的農(nóng)作物越來越多，目前廣泛用于田間生產(chǎn)的有轉(zhuǎn)基因棉花，主要用于防治棉鈴蟲。另外，轉(zhuǎn)基因水稻也可以用于二化螟、大螟和稻縱卷葉螟等害蟲的防治，同時，研究人員發(fā)現(xiàn)轉(zhuǎn)基因稻谷對印度谷螟也有不錯的防效[2]。在Bt抗蟲研究領(lǐng)域還存在幾個亟需回答的問題，一是 cry蛋白的毒性機(jī)理還未有定論，目前存在孔道假說和信號傳導(dǎo)假說兩個潛在的致死機(jī)制。另外，研究表明對印度谷螟進(jìn)行室內(nèi)篩選可以使試蟲種群獲得Bt抗性，但是相關(guān)的抗性發(fā)生機(jī)理尚未得到闡明[3]。

高通量測序技術(shù)，是核酸測序中最常用的工具，其標(biāo)志就是能一次并行對幾十萬到幾百萬條基因序列進(jìn)行測定。在高通量測序技術(shù)的迅猛發(fā)展下，生物學(xué)問題可以通過基因表達(dá)差異分析、突變分析等手段進(jìn)行解決。在昆蟲學(xué)研究領(lǐng)域，高通量測序可以用于進(jìn)化發(fā)育、生理生化以及毒理學(xué)等方面的研究。但是，印度谷螟應(yīng)用高通量測序技術(shù)的報道很少，相關(guān)基因數(shù)據(jù)有限。本研究通過對長期飼喂轉(zhuǎn)基因大米粉和常規(guī)大米粉的印度谷螟種群進(jìn)行轉(zhuǎn)錄組高通量測序（Illumina HiSeqTM 2000），組裝后分別得到了41 597和42 306條潛在基因序列，注釋其中的關(guān)鍵基因，最后比較分析差異表達(dá)基因。本研究一方面大大擴(kuò)充了該物種的基因信息，另一方面也為轉(zhuǎn)基因稻谷對印度谷螟乃至其他害蟲的抗蟲機(jī)理提供新的見解。

1 材料與方法

1.1 供試印度谷螟

本試驗(yàn)所用的印度谷螟源來源于江蘇省吳江市，后經(jīng)實(shí)驗(yàn)室多代飼養(yǎng)。其中，非脅迫品系（SS）在實(shí)驗(yàn)室條件下飼養(yǎng)20代以上（(30±0.5) ℃、相對濕度70%～80%，無光照），飼料為正常米粉（不含Bt毒素）。從非脅迫品系中隨機(jī)選擇部分試蟲用含有2%轉(zhuǎn)基因米粉的混合飼料飼養(yǎng)，篩選培養(yǎng)5個月后，視該種群為脅迫種群（RS）。

試蟲非脅迫品系（SS）飼料為：明恢63（普通稻谷）大米粉；試蟲脅迫品系（RS）飼料為：明恢63（普通稻谷）和華恢1號（Cry1Ab/Cry1Ac轉(zhuǎn)基因稻谷）大米粉混合物。上述兩種稻谷均來自華中農(nóng)業(yè)大學(xué)作物遺傳改良國家重點(diǎn)實(shí)驗(yàn)室。

1.2 RNA制備以及cDNA文庫構(gòu)建

以3齡幼蟲作為供試對象，提取RNA后送測?？俁NA提取按照RNeasy Plus Universal Mini Kit（Qiagen）試劑盒說明書進(jìn)行，用DNase I (Takara)去除基因組雜質(zhì)，測定RNA濃度，并通過1.5%瓊脂糖凝膠電泳檢測RNA完整性。然后，使用華大基因HiSeq TM 2000 (Illumina)構(gòu)建cDNA文庫并進(jìn)行測序。

1.3 生物信息學(xué)分析

篩選剔除低質(zhì)量序列后以高質(zhì)量的 reads進(jìn)行生物信息學(xué)分析。首先，Trinity組合某一重疊長度的reads以形成更長的contigs片段[4]；然后，將reads進(jìn)行末端讀數(shù)配對，確定contigs間的差異情況；最后，利用 Trinity將 contigs連接成沒有Ns的unigene序列。將組裝的unigenes與NR、Swiss-Prot、KEGG和COG數(shù)據(jù)庫進(jìn)行比對（設(shè)置臨界 E值<10–5）[5]。Unigene序列結(jié)果全部提交至 NCBI SRA（登錄號分別為 SRP060836，SRP060680）。

1.4 差異Unigenes分析

利用RPKM方法（每kb每百萬讀數(shù)）來鑒定兩個樣品（RS和SS）之間的基因表達(dá)差異[6]，同時使用FDR方法來確定P值的閾值[7]?；蝻@著差異表達(dá)的閾值為“FDR≤0.001和 log2比≥1的絕對值”。最后以GO富集和KEGG代謝通路深入分析轉(zhuǎn)錄組數(shù)據(jù)[8]。

1.5 顯著差異基因的GO分析

將DEG對應(yīng)GO數(shù)據(jù)庫，計(jì)算GO項(xiàng)基因數(shù)，分析GO功能差異，計(jì)算公式如下[10]：

N：KEGG基因數(shù)；n：N中DEG數(shù)；M：特定途徑基因數(shù)；m：M中的DEG數(shù)。

1.6 顯著差異基因的代謝通路分析

通過公共代謝路徑相關(guān)數(shù)據(jù)庫KEGG本體富集分析進(jìn)一步了解顯著差異基因的生物學(xué)功能[9]，計(jì)算公式同上[10]。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組數(shù)據(jù)De novo組裝

我們從脅迫種群（RS）和非脅迫種群（SS）中均獲得超過 51（百萬）的 reads數(shù)，且Q20＞98%，測序質(zhì)量較好，深度和覆蓋率結(jié)果如圖 1所示。之后，通過 Trinity軟件組裝 reads，共分別獲得54 648和54 647個contigs，平均長度分別均在600 bp以上。之后再進(jìn)行unigenes組裝，RS-Unigenes和SS-unigenes個數(shù)分別為41 597和42 306個，平均長度均超過1 000 bp，另外還獲得 37 246個All-unigenes（表 1）。

表1 測序產(chǎn)量統(tǒng)計(jì)

圖1 不同種群Unigene的分布情況

2.2 Unigene功能注釋與COG分類

在37 246個unigenes中有23 310個unigenes可以被注釋，其中有22 211個unigenes在NR數(shù)據(jù)庫中可以被注釋。通過NR數(shù)據(jù)庫中的E值分析，20.5%的注釋序列E值小于1.0E–100，15.5%注釋序列E值范圍在1.0E–100到1.0E–60之間，表明二者之間具有強(qiáng)同源性（表2）。序列相似性分布表明，26.1%的序列之間有較高的同源性，而另外73.9%的序列同源性介于17%～80%。另外，與印度谷螟相似度最高的物種是帝王蝶（American monarch, 67.60%），其次是家蠶（Bombyx mori, 6.5%）和赤擬谷盜（Tribolium castaneum, 3.4%）（圖 2）。

表2 注釋結(jié)果統(tǒng)計(jì)

圖2 Nr注釋的物種分布

2.3 注釋基因功能研究

10 694個序列可以注釋到 GO的 58個功能組，主要包括：生物過程，細(xì)胞組分和分子功能（圖3）。同時，細(xì)胞加工（cellular process）（6 665,62.32%）、細(xì)胞組分（cell）（5 063, 47.34%）和催化活性（catalytic activity）（5 101, 47.7%）為主。

圖3 Unigene的GO分類圖

共9 348條unigenes注釋到COG并分為25個功能類別（圖 4）。其中，最大的聚類組是“一般功能基因”（3 662, 39.17%），其次是“復(fù)制、重組和修復(fù)基因”（1 575, 16.85%）以及“翻譯、核糖體結(jié)構(gòu)和生物發(fā)生”（1 519, 16.25%）。而其它類別基因較少，如“細(xì)胞外結(jié)構(gòu)”（92, 0.98%），“RNA加工和修飾”（88, 9.41%）和“核結(jié)構(gòu)”（5,0.05%）。

圖4 COG功能注釋

KEGG注釋結(jié)果顯示，15 750條unigenes被注釋到258個KEGG通路。其中，以新陳代謝通路、肌動蛋白細(xì)胞骨架調(diào)節(jié)通路為主，另外有少數(shù)的 unigene被注釋到了生物素代謝和賴氨酸生物合成通路。

2.4 不同品系印度谷螟基因表達(dá)分析

對印度谷螟脅迫品系與非脅迫品系的差異基因進(jìn)行比較分析發(fā)現(xiàn)，有15 741條基因上調(diào)，18 725條基因下調(diào)（圖5）。進(jìn)一步利用閾值（FDR≤0.01）和log2Ratio (≥1)分析發(fā)現(xiàn)，兩個種群中差異極顯著的基因共有10 224條，其中4 520條基因上調(diào)，5 704條基因下調(diào)；除此之外，有1 245條基因在脅迫品系中表達(dá)，同時，有503條基因只在非脅迫種群中表達(dá)（圖6）。

圖5 印度谷螟不同品系差異表達(dá)的基因

圖6 印度谷螟不同品系顯著差異表達(dá)基因

2.5 印度谷螟解毒和Bt毒素受體相關(guān)基因分析

結(jié)果表明，有 87條羧酸酯酶和谷胱甘肽-S-轉(zhuǎn)移酶的相關(guān)基因在脅迫和非脅迫種群中具有表達(dá)差異，其中60.91%的基因在脅迫品系中表達(dá)量上調(diào)（表3）。

另外，我們還對不同種群試蟲的Bt毒素受體相關(guān)的基因進(jìn)行了比較分析。結(jié)果表明，有78條注釋為氨肽酶–N，堿性磷酸酶，鈣黏蛋白，糖脂類的顯著差異基因在兩個品系中被發(fā)現(xiàn)，其中有55.13%的基因在脅迫品系中上調(diào)。值得注意的是，所有注釋為 APN的基因在脅迫品系中都顯著上調(diào)。之后，我們對篩選出來的APN基因進(jìn)行了系統(tǒng)發(fā)育樹構(gòu)建，6個APN基因能夠劃到4個相關(guān)家族，與家蠶、小菜蛾、煙草天蛾等鱗翅目昆蟲的APN蛋白序列相似度較高（圖7）。

表3 印度谷螟脅迫和敏感品系中差異顯著的羧酸酯酶和谷胱甘肽-S-轉(zhuǎn)移酶基因

圖7 印度谷螟和其他鱗翅目昆蟲中APN基因系統(tǒng)發(fā)育分析

3 討論

印度谷螟危害嚴(yán)重且全球性分布，是轉(zhuǎn)基因水稻的主要靶標(biāo)害蟲之一。過去幾年，和Bt抗性機(jī)制相關(guān)的解毒代謝機(jī)制以及毒素受體被很多人研究[11-12]。然而印度谷螟對Bt毒素的抗性機(jī)理尚不清晰。本研究中，我們用二代測序技術(shù)對印度谷螟3齡敏感及經(jīng)轉(zhuǎn)基因大米粉飼喂過的脅迫品系進(jìn)行轉(zhuǎn)錄組比較分析，以研究試蟲對Bt的相關(guān)毒理學(xué)機(jī)制。通過測序，總計(jì)產(chǎn)出超過100萬nt數(shù)據(jù)，這些轉(zhuǎn)錄數(shù)據(jù)經(jīng)過過濾被組裝為Unigenes，并注釋到NR，NT，KEGG，COG和GO數(shù)據(jù)庫，其中，注釋到NR中的Unigenes最多，有22 211個，同時發(fā)現(xiàn)和印度谷螟最相似的物種是帝王蝶（67.6%），其次是家蠶（6.5%）。差異表達(dá)基因分析結(jié)果發(fā)現(xiàn)，約有46%的差異表達(dá)基因在脅迫品系中顯著上調(diào)。大量研究表明，氨肽酶-N、堿性磷酸酶、鈣黏蛋白、糖脂類可能是Bt毒素的識別受體[13-14]。在本研究，我們發(fā)現(xiàn)9個編碼APN的基因，這些基因在脅迫品系中的表達(dá)量明顯升高，另外還發(fā)現(xiàn) 55個編碼 Cadherin的基因，10個編碼ALP和4條編碼Glycolipid的基因，但這些受體基因在不同品系中的表達(dá)量無顯著差異。由此可以推測，印度谷螟對轉(zhuǎn)基因大米粉的毒理學(xué)機(jī)制和APN受體密切相關(guān)。另外，我們的另一項(xiàng)針對取食轉(zhuǎn)基因大米粉后的印度谷螟體內(nèi)抗氧化酶的活力動態(tài)進(jìn)行了深入探索，這將與本研究一起為后續(xù)的相關(guān)研究提供了堅(jiān)實(shí)的理論依據(jù)[15]。

本研究通過二代測序技術(shù)對印度谷螟轉(zhuǎn)錄組進(jìn)行了解析，同時利用生物信息分析方法對不同轉(zhuǎn)基因大米粉飼喂種群的表達(dá)差異基因進(jìn)行了比較研究，研究結(jié)果將有利于后續(xù)基因功能驗(yàn)證和印度谷螟的生物學(xué)相關(guān)研究，為儲糧害蟲防治提供了重要的理論依據(jù)和參考資料。

備注：本文的彩色圖表可從本刊官網(wǎng)（http://lyspkj.ijournal.cn/ch/index.axpx）、中國知網(wǎng)、萬方、維普、超星等數(shù)據(jù)庫下載獲取。