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

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

（南京財經大學 食品科學與工程學院，江蘇高校糧油質量安全控制及深加工重點實驗室，江蘇 南京 210023）

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

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

1 材料與方法

1.1 供試印度谷螟

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

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

1.2 RNA制備以及cDNA文庫構建

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

1.3 生物信息學分析

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

1.4 差異Unigenes分析

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

1.5 顯著差異基因的GO分析

將DEG對應GO數據庫，計算GO項基因數，分析GO功能差異，計算公式如下[10]：

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

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

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

2 結果與分析

2.1 轉錄組數據De novo組裝

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

表1 測序產量統計

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

2.2 Unigene功能注釋與COG分類

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

表2 注釋結果統計

圖2 Nr注釋的物種分布

2.3 注釋基因功能研究

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

圖3 Unigene的GO分類圖

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

圖4 COG功能注釋

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

2.4 不同品系印度谷螟基因表達分析

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

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

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

2.5 印度谷螟解毒和Bt毒素受體相關基因分析

結果表明，有 87條羧酸酯酶和谷胱甘肽-S-轉移酶的相關基因在脅迫和非脅迫種群中具有表達差異，其中60.91%的基因在脅迫品系中表達量上調（表3）。

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

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

圖7 印度谷螟和其他鱗翅目昆蟲中APN基因系統發育分析

3 討論

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

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

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