棉鈴蟲Bt抗感品系的轉(zhuǎn)錄組差異表達分析

劉運澤 張丹丹 曲愛軍 鄭方強 蕭玉濤

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

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

中圖分類號：S435.622+.3：Q786文獻標識號：A文章編號：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

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

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

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

1 材料與方法

1.1 供試昆蟲樣品制備與測序

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

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

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

1.2 轉(zhuǎn)錄本拼接

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

1.3 轉(zhuǎn)錄本定量

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

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

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

1.5 表達矩陣差異分析及注釋

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

2 結(jié)果與分析

2.1 轉(zhuǎn)錄本拼接

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

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

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

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

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

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

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

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

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

3 討論與結(jié)論

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

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

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

參 考 文 獻：

[1] 陸宴輝， 姜玉英， 劉杰， 等. 種植業(yè)結(jié)構(gòu)調(diào)整增加棉鈴蟲的災(zāi)變風險[J]. 應(yīng)用昆蟲學報， 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)棉鈴蟲技術(shù)的改進[J].植物保護，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.