劉運(yùn)澤 張丹丹 曲愛(ài)軍 鄭方強(qiáng) 蕭玉濤

摘要：棉鈴蟲(chóng)（Helicoverpa armigera）是全球性農(nóng)業(yè)害蟲(chóng)，極大地影響作物產(chǎn)量。目前各地主要通過(guò)種植Bt作物進(jìn)行防治，但近年來(lái)發(fā)現(xiàn)田間棉鈴蟲(chóng)出現(xiàn)了Bt抗性。本研究對(duì)單個(gè)棉鈴蟲(chóng)抗性品系與敏感品系進(jìn)行轉(zhuǎn)錄組測(cè)序，并使用生物信息學(xué)方法分析得到差異表達(dá)基因。結(jié)果顯示，Bt差異基因主要富集在氨肽酶、蛋白酶通路，可能對(duì)抗性產(chǎn)生影響。該研究為解釋棉鈴蟲(chóng)Bt抗性機(jī)理、挖掘新的抗性基因提供了重要依據(jù)。

關(guān)鍵詞：棉鈴蟲(chóng);Bt抗性;轉(zhuǎn)錄組測(cè)序;差異基因

中圖分類(lèi)號(hào)：S435.622+.3：Q786文獻(xiàn)標(biāo)識(shí)號(hào)：A文章編號(hào)：1001-4942（2019）07-0005-05

Abstract Cotton bollworm （Helicoverpa armigera） is a global agricultural pest that greatly affects crop yield. At present， it is mainly controlled by planting Bt crop， but Bt resistance has emerged in recent years. In this study， one laboratorial cotton bollworm resistant strains and one sensitive strain were profiled using RNA-Seq technology. With the latest database and high-cited tools， the differentially expressed genes were screened out. The results showed that the Bt-responsive genes were mainly enriched in aminopeptidase and protease pathways which might affect the resistance. This study provided an important basis for explaining the mechanism of Bt resistance in cotton bollworm and mining new resistance genes.

Keywords Helicoverpa armigera; Bt resistance; RNA-Seq; Differentially expressed genes

棉鈴蟲(chóng)（Helicoverpa armigera）是重要的多食性農(nóng)業(yè)害蟲(chóng)，可以為害棉花、玉米、小麥、蔬菜、花卉等植物，嚴(yán)重影響農(nóng)業(yè)生產(chǎn)[1]。為了有效減少棉鈴蟲(chóng)危害，從1996年至今，全球商業(yè)化種植Bt棉總面積超過(guò)9.3×108 hm2[2]。害蟲(chóng)長(zhǎng)期處于Bt毒素的高壓選擇下容易產(chǎn)生抗性，目前室內(nèi)及田間均有棉鈴蟲(chóng)Bt抗性的報(bào)道[3-5]。

近年來(lái)研究者們利用高通量轉(zhuǎn)錄組測(cè)序探索了昆蟲(chóng)差異基因、低表達(dá)基因、未知的編碼與非編碼RNA等[6，7]。對(duì)不同處理?xiàng)l件、不同發(fā)育階段或不同種類(lèi)的昆蟲(chóng)進(jìn)行轉(zhuǎn)錄組測(cè)序，可以幫助研究者更好地了解轉(zhuǎn)錄本表達(dá)豐度、基因表達(dá)差異、功能基因信息、可變剪切位點(diǎn)等關(guān)鍵信息[8-11]。Zhao等利用Illumina/Solexa平臺(tái)進(jìn)行轉(zhuǎn)錄組測(cè)序找到了棉鈴蟲(chóng)與免疫、解毒發(fā)育和代謝有關(guān)的1 022個(gè)未注釋的差異基因，進(jìn)行了免疫相關(guān)信號(hào)和調(diào)節(jié)途徑研究[12]。Wei等研究了6個(gè)不同抗性倍數(shù)的棉鈴蟲(chóng)品系的轉(zhuǎn)錄組數(shù)據(jù)，發(fā)現(xiàn)有9個(gè)胰蛋白相關(guān)基因出現(xiàn)明顯下調(diào)，并且不同品系篩選的差異基因存在較大區(qū)別，推測(cè)不同品系抗性進(jìn)化過(guò)程中可能產(chǎn)生不同的抗性機(jī)制[13]。Li 等通過(guò) Solexa 對(duì)棉鈴蟲(chóng)幼蟲(chóng)轉(zhuǎn)錄組進(jìn)行測(cè)序，獲得37 352 條 contigs，鑒定到數(shù)個(gè)可進(jìn)行 RNA 干擾的備選抗性靶基因[14]。

棉鈴蟲(chóng)等非模式昆蟲(chóng)，由于基因組缺失或注釋不完善，目前主要采用從頭拼接轉(zhuǎn)錄本的無(wú)參分析流程。本研究從生物信息分析角度出發(fā)，采用高質(zhì)量分析軟件與最新的數(shù)據(jù)庫(kù)進(jìn)行棉鈴蟲(chóng)轉(zhuǎn)錄組分析，結(jié)合功能注釋和富集分析篩選了可能與Bt抗性有關(guān)的差異表達(dá)基因，為后續(xù)的棉鈴蟲(chóng)抗性分子生物學(xué)研究提供數(shù)據(jù)支持。

1 材料與方法

1.1 供試?yán)ハx(chóng)樣品制備與測(cè)序

試驗(yàn)所用的敏感品系（susceptible strains，SS）棉鈴蟲(chóng)于1996年取自中國(guó)農(nóng)業(yè)科學(xué)院植物保護(hù)研究所河南新鄉(xiāng)基地后，一直在植物保護(hù)研究所養(yǎng)蟲(chóng)間隔離飼養(yǎng)[4]; DU抗性品系是2017年6—9月在中國(guó)農(nóng)業(yè)科學(xué)院植物保護(hù)研究所山東夏津試驗(yàn)基地利用3.6 g/L Cry1Ac毒素飼料篩選的高抗品系。飼養(yǎng)溫度為（27±2）℃，相對(duì)濕度為75%±10%，光周期為 14L∶10D。

試驗(yàn)于2018年1月進(jìn)行，首先將初孵幼蟲(chóng)接到24孔板中飼養(yǎng)，每孔4～6頭，待幼蟲(chóng)長(zhǎng)至4齡時(shí)（約孵化后第7天）移入指形管，待5齡時(shí)（約分裝后第3天）取出，準(zhǔn)備提取中腸組織。棉鈴蟲(chóng)飼養(yǎng)期間均采用正常人工飼料，其配制方法見(jiàn)梁革梅等[15]的方法。

收集SS、DU品系棉鈴蟲(chóng)的5齡幼蟲(chóng)活體，先將棉鈴蟲(chóng)在冰上放置3 min，待其不活躍后，用大頭針固定蟲(chóng)體，使用鑷子、剪刀取出中腸組織。放入4℃預(yù)冷的0.7% NaCl溶液（約10 s）洗干凈內(nèi)含物，接著用濾紙吸干水分，放入RNase-free EP管中，立刻液氮冷卻，最后-80℃保存?zhèn)溆?。每個(gè)品系3個(gè)生物學(xué)重復(fù)，每個(gè)重復(fù)15頭棉鈴蟲(chóng)。用Trizol試劑進(jìn)行總RNA提取，RNA樣品檢測(cè)合格后交由北京諾禾致源科技股份有限公司進(jìn)行測(cè)序。

1.2 轉(zhuǎn)錄本拼接

先利用FastQC、TrimGalore軟件[16]對(duì)2個(gè)品系共6個(gè)樣本的原始數(shù)據(jù)進(jìn)行質(zhì)量控制，保證不含有低質(zhì)量堿基及測(cè)序接頭，然后利用Trinity軟件[17]對(duì)6個(gè)樣本混合拼接，對(duì)拼接結(jié)果進(jìn)行質(zhì)量控制，為保證準(zhǔn)確性同時(shí)采用了4種方式：計(jì)算N50值;回貼比對(duì);利用BLAST[18]、HMMER[19]結(jié)合BUSCO[20]數(shù)據(jù)庫(kù)分析組裝結(jié)果中的同源基因;使用Diamond軟件[21]比對(duì)非冗余蛋白數(shù)據(jù)庫(kù)（Nr）。

1.3 轉(zhuǎn)錄本定量

采用不基于比對(duì)的Salmon軟件[22]將2個(gè)品系的6個(gè)樣本分別定量然后組合成一個(gè)表達(dá)矩陣。

1.4 組內(nèi)樣品重復(fù)性評(píng)估和組間聚類(lèi)

對(duì)利用Salmon得到的表達(dá)矩陣基于R語(yǔ)言計(jì)算每組生物學(xué)重復(fù)之間的Pearson相關(guān)系數(shù)，并繪制熱圖比較組間與組內(nèi)的生物學(xué)重復(fù)相關(guān)性。

1.5 表達(dá)矩陣差異分析及注釋

利用edgeR[23]分析表達(dá)矩陣的原始統(tǒng)計(jì)值，利用calcNormFctors方法進(jìn)行組間標(biāo)準(zhǔn)化，得到差異表達(dá)矩陣，接著利用其中的log2（fold change）列（以下簡(jiǎn)稱(chēng)logFC），將火山圖代碼中的logFC閾值設(shè)為mean（abs（logFC））+2×sd（abs（logFC）），即logFC絕對(duì)值的均值加兩倍絕對(duì)值標(biāo)準(zhǔn)差，設(shè)置pvalue和qvalue閾值為0.05。計(jì)算DU-SS比較組的上調(diào)、下調(diào)差異基因。之后利用clusterProfiler的R包結(jié)合最新的GO、KEGG數(shù)據(jù)庫(kù)進(jìn)行功能注釋和富集分析，結(jié)果中篩選出pvalue和qvalue小于0.05的GO條目與KEGG 通路進(jìn)行統(tǒng)計(jì)和分析。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄本拼接

SS與DU品系的測(cè)序數(shù)據(jù)結(jié)果均得到2 000×104以上的reads，平均72×108 bp的轉(zhuǎn)錄組數(shù)據(jù)，過(guò)濾后堿基質(zhì)量均達(dá)到Q30以上。對(duì)利用Trinity拼接得到的166 594條轉(zhuǎn)錄本進(jìn)行質(zhì)量檢測(cè)，N50為1 036 bp。將SS與DU品系的測(cè)序回貼比對(duì)到拼接好的轉(zhuǎn)錄本，得到SS比對(duì)率為98.28%，DU比對(duì)率為98.35%，說(shuō)明拼接得到的轉(zhuǎn)錄本中低表達(dá)、低質(zhì)量或重復(fù)reads很少，拼接效果很好。BUSCO檢測(cè)到昆蟲(chóng)保守基因的完整度為95%，一般認(rèn)為這個(gè)值大于80%表示轉(zhuǎn)錄本拼接質(zhì)量是可靠的。Nr數(shù)據(jù)庫(kù)比對(duì)結(jié)果中，棉鈴蟲(chóng)有47%的轉(zhuǎn)錄本序列比對(duì)到棉鈴蟲(chóng)自身，其次與煙草夜蛾相似性達(dá)到了21%，另外還與家蠶、斜紋夜蛾、二化螟等鱗翅目昆蟲(chóng)序列較為相似（圖1）。

2.2 組內(nèi)樣品重復(fù)性評(píng)估和組間聚類(lèi)

對(duì)得到的表達(dá)矩陣，計(jì)算每組3個(gè)生物學(xué)重復(fù)之間的相關(guān)性。一般而言，樣本間重復(fù)的Pearson相關(guān)系數(shù)大于0.8即為存在相關(guān)性。結(jié)果顯示SS與DU品系棉鈴蟲(chóng)各自的生物學(xué)重復(fù)之間Pearson系數(shù)均在0.9附近，表示制備的生物學(xué)重復(fù)樣本合格（圖2A、B）。另外組間聚類(lèi)顯示SS與DU樣本各自聚在一起，且表達(dá)差異明顯，說(shuō)明測(cè)序結(jié)果符合試驗(yàn)設(shè)計(jì)預(yù)期（圖2C）。

2.3 轉(zhuǎn)錄本差異分析

利用edgeR對(duì)Salmon表達(dá)矩陣進(jìn)行差異分析，共得到1 156個(gè)差異基因，其中上調(diào)表達(dá)有593個(gè)，下調(diào)表達(dá)有563個(gè)（圖3）。火山圖中顯示的上調(diào)及下調(diào)“離群點(diǎn)”表示DU-SS比較組中差異顯著性更強(qiáng)的基因，這些差異極顯著的基因更有可能導(dǎo)致抗性的產(chǎn)生。

2.4 差異基因功能注釋與富集分析

DU品系與SS品系得到的1 156個(gè)差異表達(dá)基因中，有397個(gè)可以比對(duì)到Nr庫(kù)并對(duì)應(yīng)到RefSeq數(shù)據(jù)庫(kù)。對(duì)其進(jìn)行GO分析，發(fā)現(xiàn)了66個(gè)條目顯著性富集。其中GO 生物學(xué)過(guò)程包括36個(gè)條目，細(xì)胞組分（cellular component，CC）包括10個(gè)條目，分子功能（molecular function，MF）包括20個(gè)條目。按照pvalue值從小到大挑選前20個(gè)條目進(jìn)行可視化（圖4A）。在分子功能這一類(lèi)中，有5個(gè)基因（HaOG204148、 HaOG204182、LOC110378622、HaOG212290、HaOG216726）與氨肽酶活性有關(guān)（GO：0004177），氨肽酶如果與水解后的Bt毒素單體結(jié)合會(huì)導(dǎo)致中腸穿孔，因此氨肽酶活性可以作為判斷抗性的指標(biāo)，也說(shuō)明本研究篩選得到的差異表達(dá)基因具有一定的可靠性。

另有12個(gè)基因（HaOG200516、HaOG200394、HaOG204148、HaOG204182、HaOG200489、HaOG200480、HaOG203137、HaOG215606、LOC110378622、HaOG212290、HaOG200476、HaOG216726）會(huì)影響蛋白酶活性（GO：0008233），蛋白酶的缺失會(huì)使Bt前毒素活化程度降低，導(dǎo)致棉鈴蟲(chóng)對(duì)Bt毒素產(chǎn)生抗性。此外，有2個(gè)基因（HaOG216480、HaOG210130）響應(yīng)抗菌免疫過(guò)程（GO：0002812），推測(cè)免疫反應(yīng)與Bt抗性存在一定的聯(lián)系。

KEGG結(jié)果同樣選取pvalue前20個(gè)通路進(jìn)行可視化（圖4B）。其中有10個(gè)基因（[JP3]HaOG212290、HaOG205566、HaOG215606、HaOG214145、 HaOG209294、HaOG208605、LOC110378622、HaOG200470、HaOG212255）在肽酶通路（ko01002）中發(fā)揮作用，支持了GO中因氨肽酶活性變化而產(chǎn)生抗性的推測(cè)。

3 討論與結(jié)論

棉鈴蟲(chóng)長(zhǎng)期處于Bt毒蛋白的高壓選擇下，容易產(chǎn)生抗性，Bt抗性機(jī)理研究一直是防治棉鈴蟲(chóng)的關(guān)鍵一步。中腸是昆蟲(chóng)食物消化與吸收的主要場(chǎng)所，也是Bt毒素發(fā)揮作用的重要組織。本研究利用以3.6 g/L的Bt毒飼料培養(yǎng)出的抗性品系，通過(guò)中腸轉(zhuǎn)錄組測(cè)序分析得到該抗性品系與敏感品系的差異基因。

本研究篩選到的棉鈴蟲(chóng)抗感品系差異基因主要與氨肽酶、蛋白酶活性有關(guān)，參與了肽酶通路。前人研究證實(shí)Cry1Ac能與氨肽酶受體的N-乙酰氨基半乳糖（GalNAc）基團(tuán)結(jié)合[24]，并且氨肽酶的突變或表達(dá)量降低都會(huì)增加棉鈴蟲(chóng)對(duì)Cry1Ac毒素的抗性[25]。中腸蛋白酶與Cry前體毒素活化關(guān)系密切。蛋白酶活性降低導(dǎo)致前毒素活化程度不足，不能產(chǎn)生足夠的活化毒素，因此會(huì)產(chǎn)生抗性。已有研究發(fā)現(xiàn)煙芽夜蛾（Heliothis virescens）和歐洲玉米螟（Ostrinia nubilalis）的抗性品系中參與前毒素活化的中腸蛋白酶的活性下降[26，27]。棉鈴蟲(chóng)中腸缺乏絲氨酸蛋白酶也會(huì)導(dǎo)致前毒素活化不足，從而產(chǎn)生Bt抗性[28]。另外有2個(gè)篩選到的基因可以影響抗菌免疫反應(yīng)，前人在研究地中海粉螟（Ephestia kuehniella）時(shí)發(fā)現(xiàn)蟲(chóng)體內(nèi)免疫反應(yīng)會(huì)產(chǎn)生凝聚物，它可以與Cry1Ac結(jié)合，導(dǎo)致Cry1Ac與中腸受體蛋白結(jié)合減少[29]，說(shuō)明長(zhǎng)期處于高濃度Bt毒素條件下的棉鈴蟲(chóng)免疫反應(yīng)較強(qiáng)。以上分析的差異基因未曾報(bào)道過(guò)，它們?cè)诿掴徬x(chóng)抗性過(guò)程中的功能需要結(jié)合分子試驗(yàn)進(jìn)一步驗(yàn)證。

本研究結(jié)果為進(jìn)一步篩選與克隆關(guān)鍵Bt抗性基因提供了依據(jù)，也為探究棉鈴蟲(chóng)Bt抗性機(jī)制奠定了基礎(chǔ)。

參 考 文 獻(xiàn)：

[1] 陸宴輝， 姜玉英， 劉杰， 等. 種植業(yè)結(jié)構(gòu)調(diào)整增加棉鈴蟲(chóng)的災(zāi)變風(fēng)險(xiǎn)[J]. 應(yīng)用昆蟲(chóng)學(xué)報(bào)， 2018， 55（1）： 19-24.

[2] Abbas M S T. Genetically engineered （modified） crops （Bacillus thuringiensis crops） and the world controversy on their safety[J]. Egyptian Journal of Biological Pest Control， 2018， 28（1）： 28-52.

[3] Tabashnik B E， Wu K M， Wu Y D. Early detection of field-evolved resistance to Bt cotton in China： Cotton bollworm and pink bollworm[J]. Journal of Invertebrate Pathology， 2012， 110（3）： 301-306.

[4] Guo Y Y， Li K K， Wu K M， et al. Changes of inheritance mode and fitness in Helicoverpa armigera（Hübner） （Lepidoptera∶Noctuidae） along with its resistance evolution to Cry1Ac toxin[J]. Journal of Invertebrate Pathology， 2007， 97（2）： 142-149.

[5] Downes S， Mahon R. Successes and challenges of managing resistance in Helicoverpa armigera to Bt cotton in Australia[J]. GM Crops and Food， 2012， 3（3）： 228-234.

[6] Qazi S S， Mittapalli O， Beliveau C， et al. Identification of odor-processing genes in the emerald ash borer，Agrilus planipennis[J]. PLoS One， 2013， 8（2）： e56555.

[7] Legeai F， Derrien T. Identification of long non-coding RNAs in insects genomes[J]. Current Opinion in Insect Science， 2015， 7： 37-44.

[8] Li J， Zhu L， Hull J J， et al. Transcriptome analysis reveals a comprehensive insect resistance response mechanism in cotton to infestation by the phloem feeding insect Bemisia tabaci（whitefly）[J]. Plant Biotechnology Journal， 2016， 14（10）： 1956-1975.

[9] Chang J C， Ramasamy S. Transcriptome analysis in the beet webworm， Spoladea recurvalis（Lepidoptera∶Crambidae）[J]. Insect Science， 2018， 25（1）： 33-44.

[10]Ross H A， Buckley T R， Dennis A B， et al. De Novo Transcriptome analysis of the common new Zealand stick insect Clitarchus hookeri（Phasmatodea） reveals genes involved in olfaction， digestion and sexual reproduction[J]. PLoS One， 2016， 11（6）： e0157783.

[11]Katsuma S， Kawamoto M， Shoji K， et al. Transcriptome profiling reveals infection strategy of an insect maculavirus[J]. DNA Research， 2018， 25（3）： 277-286.

[12]Zhao Z， Wu G， Wang J， et al. Next-generation sequencing-based transcriptome analysis of Helicoverpa armigera larvae immune-primed with Photorhabdus luminescens TT01[J]. PLoS One， 2013， 8（11）： e80146.

[13]Wei J， An S， Yang S， et al. Transcriptomice responses to different Cry1Ac selection stresses in Helicoverpa armigera[J]. Frontiers in Physiology， 2018， 9： 1653-1670.

[14]Li J， Li X， Chen Y， et al. Solexa sequencing based transcriptome analysis of Helicoverpa armigera larvae[J]. Molecular Biology Reports， 2012， 39（12）： 11051-11059.

[15]梁革梅，譚維嘉，郭予元. 人工飼養(yǎng)棉鈴蟲(chóng)技術(shù)的改進(jìn)[J].植物保護(hù)，1999（2）：16-18.

[16]Marcel M. Cutadapt removes adapter sequences from high-throughput sequencing reads[J]. EMBnet. Journal， 2011， 17（1）： 10-12.

[17]Grabherr， Manfred G， Brian J，et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome [J]. Nature Biotechnology， 2011， 29（7）： 644-652.

[18]Madden T L， McGinnis S. BLAST： at the core of a powerful and diverse set of sequence analysis tools[J]. Nucleic Acids Research， 2004， 32： 20-25.

[19]Johnson L S， Eddy R， Portugaly E. Hidden Markov ?model speed heuristic and iterative HMM search procedure[J]. BMC Bioinformatics， 2010， 11： 431-439.

[20]Simo F A，Waterhouse R M，Ioannidis P，et al.BUSCO：assessing genome assembly and annotation completeness with single-copy orthologs[J]. Bioinformatics， 2015， 31（19）： 3210-3212.

[21]Buchfink B， Xie C， Huson ?H. Fast and sensitive protein alignment using DIAMOND[J]. Nature Methods， 2014， 12（1）： 59-60.

[22]Patro R， Duggal G， Love M I， et al. Salmon provides fast and bias-aware quantification of transcript expression[J]. Nature Methods， 2017， 14（4）： 417-419.

[23]Robinson M D， McCarthy D J， Smyth G K. edgeR：a Bioconductor package for differential expression analysis of digital gene expression data[J]. Bioinformatics， 2009， 26（1）： 139-140.

[24]Sarkar A， Hess D， Mondal H A， et al. Homodimeric alkaline phosphatase located at Helicoverpa armigera Midgut， a putat ive receptor of Cry1Ac contains α-GaINAc in term inal glycan structure as interactive Epitope[J]. Journal of Proteome Research， 2009， 8（4）： 1838-1848.

[25]Liang G， Wang G， Wu K， et al. Mutation of an aminopeptidase N gene is associated with Helicoverpa armigera resistance to Bacillus thuringiensis Cry1Ac toxin[J]. Insect Biochemistry and Molecular Biology， 2009， 39（7）： 421-429.

[26]Oppert B，Higgins R A，Huang F，et al.Comparative analysis of proteinase activities of Bacillus thuringiensis-resistant and -susceptible Ostrinia nubilalis（Lepidoptera∶Crambidae）[J]. Insect Biochemistry and Molecular Biology， 2004， 34（8）：753-762.

[27]Huang F， Oppert B， Higgins R A， et al. Susceptibility of Dipel-resistant and -susceptible Ostrinia nubilalis（Lepidoptera∶Crambidae） to individual Bacillus thuringiensis protoxins[J]. Journal of Economic Entomology， 2005， 98（4）： 1333-1340.

[28]Liu C X， Xiao Y T， Li X C， et al. Cis-mediated down-regulation of a trypsin gene associated with Bt resistance in cotton bollworm[J]. Scientific Reports， 2014， 4（1）： 7219-7225.

[29]Mahbubur R M， Roberts H S， Schmidt O. Tolerance to Bacillus thuringiensis endotoxin in immune-suppressed larvae of the flour moth Ephestia kuehniella[J]. Journal of Invertebrate Pathology， 2007， 96（2）：125-132.