瞿貴軍，林毅

(華僑大學 化工學院，福建 廈門 361021)

桃蚜(Myzuspersicae)的寄主植物多達50多科400余種，桃蚜主要以果樹、蔬菜、煙草、花卉等數(shù)百種植物為食，是世界性的重要經濟害蟲之一[1]，也是世界傳播病毒昆蟲之首[2].桃蚜的次級共生菌主要是保護其宿主免受病原體和天敵的侵害，增強對殺蟲劑的抗性、新陳代謝和生物合成等[3].桃蚜的生命周期復雜，目前主要通過化學殺蟲劑殺桃蚜，但桃蚜對使用過的殺蟲劑都能產生不同程度的耐藥性，具有很強的遺傳變異性和適應性[4].目前，桃蚜在抗性克隆中通過細胞色素P450(P450)和UDP葡萄糖醛酸轉移酶(UGT)高度上調表達[5].環(huán)境友好型生物農藥Cry蛋白已被國際抗性行動委員會(IRAC)分類方案收錄，作為殺蟲劑活性成分[6]，在生物農藥研發(fā)中，Cry41-related蛋白可通過增強組織蛋白酶B活性殺死桃蚜[7]，而Cry1Cb2蛋白也有望成為一種有前途的抗蚜蛋白[8].

近年來，在基因挖掘和調控機制等方面的研究中，加權基因共表達網絡分析(WGCNA)[9]、差異表達基因分析(DEGs)[10]、GO/KEGG富集分析、蛋白質互作網絡分析(PPI)等綜合生物信息技術分析取得了許多研究成果.He等[11]通過WGCNA揭示昆蟲真菌病原體的各種不利脅迫適應性響應的基因和信號通路.Li等[12]通過差異表達基因分析，為評估傳毒和非傳毒桃蚜對植物病毒的反應機制提供依據(jù).Chin等[13]使用Cytoscape軟件中的CytoHubba插件，根據(jù)網絡特征對網絡中的節(jié)點進行排名.基于此，本文對GSE18659，GSE37310等2個GEO數(shù)據(jù)(16個樣本)[14-15]采用WGCNA，DEGs，PPI，同源建模及Z-dock分子對接等生物信息技術進行桃蚜關鍵抗性基因的挖掘及抗蚜Cry蛋白預測.

1 材料與方法

1.1 基因芯片的重注釋

在NCBI網站(https:∥www.ncbi.nlm.nih.gov)下載桃蚜的基因芯片GPL9470探針序列文件和G0061.0版基因組注釋文件，采用序列比對工具SeqMap(http:∥soft.sangerbox.com)和基因注釋工具，通過序列比對的方式，將基因組重新比對到探針序列上，從而獲得最新的基因編號和注釋.

1.2 關鍵抗性基因的挖掘

首先，采用加權基因共表達網絡分析，在Sangerbox網站(http:∥sangerbox.com/Too)上對桃蚜抗性研究的GSE18659，GSE37310等2個GEO數(shù)據(jù)(16個樣本)進行在線分析，離群樣本閾值設定為0.00，共表達網絡類型選擇unsign，hub相關性閾值設定為0.9，網絡權重閾值設定為0.00，模塊最小基因數(shù)設定為20.00，模塊合并閾值設定為0.25，敏感性(1～4)設定為2.00，軟閥值設定為0.00，獲得關鍵抗性基因文件和基因互作文件.

然后，采用Cytoscape軟件的Cytohubba插件，對已用R語言篩選出weight值大于0.48的基因互作文件進行打分，通過MCC算法和Degree算法，得到兩種算法各前10名基因探針對應的基因.

最后，通過GEO數(shù)據(jù)庫自帶基于limma包R語言的GEO2R工具，對GSE18659進行差異表達基因分析，根據(jù)校正后的P值Padj<0.05，差異倍數(shù)|log FC| >1，篩選出差異表達的關鍵抗性基因.在DAVID網站(https:∥david.ncifcrf.gov)選擇桃蚜研究的模式生物豌豆蚜Buchnera aphidicola str.APS (Acyrthosiphonpisum)生物體進行在線分析，獲得GO/KEGG富集分析數(shù)據(jù).

1.3 抗性調控網絡的構建

將挖掘到的關鍵抗性基因通過String數(shù)據(jù)庫(https:∥string-db.org)構建抗性調控網絡.將模式生物豌豆蚜(Acyrthosiphonpisum)的蛋白質互作數(shù)據(jù)(String數(shù)據(jù)庫)與桃蚜數(shù)據(jù)(Uniprot數(shù)據(jù)庫(https:∥www.uniprot.org))進行合并，在豌豆蚜中去掉與桃蚜名稱相同的蛋白，再與WGCNA和DEGs獲得的桃蚜關鍵抗性基因進行名稱比對，根據(jù)與關鍵基因描述名稱是一致的原則進行篩選，獲得蛋白質互作網絡分析所需的蛋白名.

在String數(shù)據(jù)庫中進行蛋白質互作分析時，分別選擇桃蚜共生菌Buchnera aphidicola str.USDA(Myzuspersicae)生物體和豌豆蚜(Acyrthosiphonpisum)生物體，通過桃蚜和豌豆蚜的相同蛋白名互換基因名實現(xiàn)兩個調控網絡串聯(lián)，合并繪制桃蚜抗性調控網絡.同時，通過NCBI網站提供的blastx工具，將前期WGCNA和DEGs得到的關鍵抗性基因序列轉化為蛋白序列，在String數(shù)據(jù)庫選擇桃蚜共生菌Buchnera aphidicola str.USDA(Myzuspersicae)生物體進行蛋白質互作分析.將經Multiple proteins和Multiple sequences檢索策略得到的蛋白互作數(shù)據(jù)合并，得到string_interactions.tsv文件.最后，采用Cytoscape軟件打開該數(shù)據(jù)重新繪制，獲得桃蚜抗性調控網絡圖；采用基迪奧生物云工具(https:∥www.omicshare.com)繪制抗性調控權重網絡圖，連通性方法選擇硬閥值，權重范圍選擇默認值(0.5～1.0).

1.4 關鍵靶標蛋白庫的構建

首先，在桃蚜抗性調控網絡中篩選出DEGs中上調表達的關鍵基因.

然后，通過Cytoscape軟件的ClueGO，CluePedia插件，選擇模式生物豌豆蚜生物體對桃蚜關鍵抗性基因進行分析，得到值得關注的蛋白.

最后，整合蛋白互作數(shù)據(jù)的分析結論，繪制Cry蛋白殺桃蚜的靶標蛋白調控網絡.采用NCBI網站的Blastp工具，在UniProt網站輸入靶標蛋白名或登錄號，得到同源序列.利用Swiss-MODEL網站(https:∥swissmodel.expasy.org)進行同源建模，選擇“Experimental Structures”的Range覆蓋率最高或“Homology model”的QMEANDisCo分數(shù)最高，確定靶標蛋白庫的分子對接實驗模型.

1.5 Cry抗蚜蛋白的預測

在Molprobity網站(http:∥molprobity.biochem.duke.edu/index.php)對同源建模的3D模型結果進行評估，在Z-dock網站(https:∥zdock.umassmed.edu)將評估合格的靶標蛋白與Cry蛋白的3D結構進行在線分子對接，所有參數(shù)設置為默認值，將實驗結果文件導入PyMol軟件進行觀察，并與參考組進行對比.

2 實驗結果與分析

2.1 桃蚜關鍵抗性基因的挖掘

通過Cytohubba插件得到MCC算法和Degree算法的前10個關鍵抗性基因探針，結果如表1所示.由表1可知：MCC算法的第3位探針與Degree算法的第5位探針相同；MCC算法的第7位探針與Degree算法的第2位探針相同；MCC算法的第10位探針與Degree算法的第10位探針相同，且均為未注釋的基因及蛋白；Degree算法的第8位探針是腺苷三磷酸(ATP)依賴性蛋白The DEAD-box RNA helicase.

表1 MCC算法和Degree算法的前10個關鍵抗性基因探針Tab.1 Top 10 key resistance gene probes of MCC algorithm and Degree algorithm

通過加權基因共表達網絡分析，可得26個基因共表達模塊，將基因探針通過Gene symbol去重后，得到2 426個樞紐基因.抗性蛋白在模塊中的分布情況，如表2所示.

由表2可知：細胞色素P450、UDP葡萄糖醛酸轉移酶、羧酸酯酶(CarE)和谷胱甘肽S-轉移酶(GST)主要集中于plum1，blue，cyan，brown等4個基因共表達模塊.

表2 抗性蛋白在模塊中的分布情況Tab.2 Distribution of resistance proteins in modules

模塊特征向量聚類分析圖，如圖1所示.由圖1及相關分析可知：plum1與brown模塊呈正相關，與blue，brown等模塊均呈負相關，blue模塊與cyan模塊呈高度正相關，P450與CarE在某些情況下可能存在協(xié)同表達，UGT與GST則存在協(xié)同表達；在cyan模塊中發(fā)現(xiàn)了Cathepsin B；在blue模塊中發(fā)現(xiàn)了The DEAD-box RNA helicase 52，The DEAD-box RNA helicase，activator of 90 kDa heat shock protein ATPase homolog 1，Hsp60 protein，Heat shock protein cognate 4，ATP-dependent 6-phosphofructokinase (PFKA)和6-phosphofructo-2-kinase等；在brown模塊中發(fā)現(xiàn)了ATP-binding cassette sub-family，heat shock 70 kDa protein 4和heat shock 70 kDa protein 3.

圖1 模塊特征向量聚類分析圖Fig.1 Cluster analysis diagram of module feature vectors

在差異表達基因分析中，共篩選出2 263個差異表達基因探針；在114個上調表達基因中發(fā)現(xiàn)了P450，CarE和GST，Hsp60 protein，Hsp83 protein(與分子伴侶Hsp90同源)，The DEAD-box RNA helicase 52 (DDX52)，Cathepsin B，6-phosphofructo-2-kinase.對DEGs上調表達的關鍵基因與WGCNA的關鍵樞紐基因取交集，以“Gene description”為基準去重，共獲得558個關鍵抗性基因.

2.2 關鍵抗性基因的GO/KEGG富集分析

在DAVID網站選擇模式豌豆蚜共生菌Buchnera aphidicola str.APS (Acyrthosiphonpisum)生物體進行在線分析，獲得GO/KEGG富集分析的數(shù)據(jù).與桃蚜共生菌相關的關鍵抗性基因GO富集分析，如表3所示.表3中：n為基因數(shù)；η為基因占比.由表3可知：成功轉化了58個關鍵抗性基因，在生物學過程中，關鍵抗性基因參與細胞氧化還原穩(wěn)態(tài)和糖酵解過程；在細胞定位中，關鍵抗性基因富集于細胞質；從分子功能上看，關鍵抗性基因具有使ATP結合、氨基?；鵷RNA編輯活性和核酸結合的作用.由于這些關鍵抗性基因與桃蚜抗性相關，從基因數(shù)和基因占比來看，其主要富集于細胞質和ATP結合.

表3 與桃蚜共生菌相關的關鍵抗性基因GO富集分析Tab.3 GO enrichment analysis of key resistance genes associated with symbiotic bacteria of Myzus persicae

與桃蚜共生菌相關的關鍵抗性基因KEGG富集分析，如表4所示.由表4可知：這些關鍵抗性基因主要富集于氨酰tRNA生物合成、碳代謝、不同環(huán)境中的微生物代謝、糖酵解/糖異生、RNA降解、谷胱甘肽代謝、磷酸戊糖途徑、檸檬酸循環(huán)(TCA循環(huán))等通路；pfka，eno和pyka等3個基因參與碳代謝、微生物代謝和糖酵解/糖異生等3個通路，而pfka和eno基因也參與RNA降解通路.

表4 與桃蚜共生菌相關的關鍵抗性基因KEGG富集分析Tab.4 KEGG enrichment analysis of key resistance genes associated with symbiotic bacteria of Myzus persicae

2.3 桃蚜抗性調控網絡的構建

圖2 桃蚜抗性調控網絡Fig.2 Resistance control network of Myzus persicae

桃蚜抗性調控網絡，如圖2所示.由圖2可知：桃蚜抗性調控網絡共有154個基因；最內圈黃色為上調表達的acee，dead，dnaj，acpp，actb-348，metk，sped，trxb，ydhd，sda，pyka等11個基因，這11個基因對應的桃蚜蛋白名和豌豆蚜蛋白名是一致的；中間圈為上調表達的基因；最外圈為非上調表達的基因.

2.4 關鍵靶標蛋白庫的構建

上調表達的桃蚜關鍵抗性蛋白的PPI網絡，如圖3所示.由圖3可知：在65個上調表達的關鍵抗性基因中，dead直接作用于dnak.

Cry蛋白殺桃蚜的關鍵靶標蛋白，如表5所示.

圖3 上調表達的桃蚜關鍵抗性蛋白互作網絡Fig.3 PPI network of upregulated key resistance proteins of Myzus persicae

表5 Cry蛋白殺桃蚜的關鍵靶標蛋白Tab.5 Key target proteins of Cry proteins against Myzus persicae

圖4 Cry蛋白殺桃蚜的靶標蛋白網絡Fig.4 Target proteins network of Cry proteins against Myzus persicae

對實驗用到的acpp，eno，pfka，dnak，dead，100164879，acypi006185-pa，catb-348，hsp60，100169264，acypi007750-pa等11個靶標基因繪制Cry蛋白殺桃蚜的靶標蛋白網絡，如圖4所示.由圖4可知：eno和acpp直接與pfka互作，dead能直接作用于dnak.

ClueGo術語/途徑網絡基因分布視圖，如圖5所示.由圖5可知：豌豆蚜LOC100169264是Glutathione S-transferase D7蛋白的基因，它和Cathepsin B一同參與相關蛋白降解途徑溶酶體，而Cathepsin B蛋白還參與了蛋白降解途徑自噬，Glutathione S-transferase蛋白還參與了氧化磷酸化和吞噬體等途徑；豌豆蚜LOC100160371是UDP-glucuronosyltransferase蛋白的基因，它和Glutathione S-transferase蛋白一同參與了吞噬體途徑；豌豆蚜LOC100160659是PFKA的基因，它參與了果糖和甘露糖代謝、半乳糖代謝、糖酵解/糖異生、戊糖磷酸途徑，和tRNA(鳥嘌呤(37)-N1)-甲基轉移酶一同參與了RNA降解；豌豆蚜LOC100162014是Protein CLP1 homolog的基因，它和pfka一同參與了糖酵解/糖異生途徑.

圖5 ClueGo術語/途徑網絡基因分布視圖Fig.5 Genes distribution view of ClueGo terms/pathways network

2.5 基于分子對接篩選潛在的抗蚜Cry蛋白

79個Cry蛋白與Acyl carrier protein (ACP)在Z-dock網站進行分子對接實驗的結果表明，Cry1Aa，Cry4Ba，Cry1Ba，Cry5Ba，Cry1Ka，Cry1Bc，Cry1Be，Cry1Bb，Cry1Af，Cry15Aa等10個Cry蛋白能與其發(fā)生連接.然后，將這10個Cry蛋白與Actin，APN，PFKA，Cathepsin B，Dead-CshA，DNAK等7個靶標蛋白在Z-dock網站進行分子對接實驗，參考抗蚜Cry1Cb2蛋白與這7個靶標蛋白的參考組對接規(guī)則，用PyMol軟件查看實驗組的對接結果，均能發(fā)生連接.

3 討論與分析

桃蚜抗性機制可以分為靶標敏感性降低和抗性代謝能力增強等兩大類，通過其羧酸酯酶的過量表達、乙酰膽堿酯酶的突變、電壓門控鈉通道突變、GABA受體亞基基因的復制、細胞色素P450CYP6CY3的過表達、煙堿乙酰膽堿受體(NAChR)的突變和ABC轉運蛋白等發(fā)揮蚜蟲的抗性機制[16-22]，文中對已報道的桃蚜關鍵抗性蛋白UGT 和GST預測存在協(xié)同表達，P450和CarE在某些通路下協(xié)同表達，DDX52，Cathepsin B，PFKA，Hsp60 protein等與UGT 和GST可能存在協(xié)同表達，數(shù)據(jù)顯示P450，CarE，GST，HSP60，DDX52，Cathepsin B等均以上調表達參與抗性.GO富集分析表明，關鍵抗性基因富集在細胞氧化還原穩(wěn)態(tài)和糖酵解過程，從基因數(shù)和基因占比來看，其主要富集于細胞質和ATP結合，說明桃蚜共生菌參與抗性主要是在細胞質，通過與ATP結合相關的基因調控、代謝通路或信號轉導等來實現(xiàn).KEGG富集分析表明，eno，pfka等基因參與RNA降解通路，而eno和acpp直接與pfka互作，dead能直接作用于dnak.

除已報道的UGT 和P450等桃蚜關鍵抗性蛋白，在基于String數(shù)據(jù)庫預測構建的桃蚜抗性調控網絡中，烯醇化酶(Enolase)，PFKA，CLP，DNAK，CshA，Pescadillo，Autophagy protein和ACP等蛋白也值得重點關注.烯醇化酶是由蚜蟲ERV畸胎細胞釋放，介導參與病原體和腫瘤細胞侵入組織及逃避宿主免疫反應的酶的激活[23].Ae-ENO是由胚胎膜分離的細胞產生的，這些細胞和毒液一起在宿主蚜蟲的調節(jié)和營養(yǎng)開發(fā)中發(fā)揮重要作用.PFKA是糖酵解最重要的調節(jié)酶，它和LOC100162014一同參與了糖酵解/糖異生途徑，LOC100162014是模式生物豌豆蚜幾丁質酶樣蛋白(CLP)的基因[24].此外，Cry-related蛋白能夠與桃蚜共生菌PFKA結合.DNAK是豌豆蚜細胞內共生熱休克蛋白同系物[25]，DNAK伴侶蛋白也是蚜蟲蜜露的蛋白質之一[26]，它可能在植物-蚜蟲相互作用中充當介質.CshA是一種二聚體 DEAD-box解旋酶，可與核糖核酸酶合作進行mRNA 轉換[27].上調表達的acypi006185-pa是模式生物豌豆蚜的Pescadillo homolog基因，Pescadillo是一種在惡性細胞中異常表達的細胞周期調節(jié)蛋白，在細胞增殖中起著至關重要的作用，并且可能是致癌轉化和腫瘤進展所必需的[28]，研究發(fā)現(xiàn)DEAD-box解旋酶DDX27可以與Pes1特異性相互作用[29].模式生物豌豆蚜基因10016879在String數(shù)據(jù)庫注釋是自噬蛋白，參與自噬囊泡的形成，比較轉錄組分析高粱甘蔗桃蚜抗性的遺傳機制發(fā)現(xiàn)在SCA感染后抗性基因型中自噬的基因上調表達[30]，另有研究發(fā)現(xiàn)，選擇性自噬通過NBR1介導的病毒衣殼蛋白和顆粒靶向抑制花椰菜花葉病毒侵染，病毒觸發(fā)的自噬很大程度上獨立于 NBR1 的方式，防止受感染植物的廣泛衰老和組織死亡，這種生存功能顯著延長了病毒生存的時間，從而增加了蚜蟲載體和 CaMV 傳播獲得病毒顆粒的機會[31].Autophagy protein極有可能參與了桃蚜的抗性作用.acpp是?；d體蛋白ACP的基因，ACP是一類具有保守的絲氨酸殘基的小分子量的酸性蛋白質，它也可作為抗菌劑的靶點[32]，參與桃蚜的代謝途徑和次級代謝產物的生物合成途徑，也有研究發(fā)現(xiàn)delta914:0-酰基載體蛋白脂肪酸去飽和酶基因的表達對于在抗蟲天竺葵(Pelargoniumxhortorum)中發(fā)現(xiàn)的omega 5漆樹酸是必需的[33]，參考受到廣泛認可的Bravo穿孔模型[34]，Cry蛋白由于高pH值和還原條件而溶解于中腸腔，并被腸道蛋白酶激活，產生毒素片段，但目前多數(shù)Cry蛋白殺桃蚜效果并不理想，這也可能與上調表達的桃蚜acpp基因有關系.ACP和PFKA存在蛋白質互作，而PFKA參與了碳代謝、微生物代謝和糖酵解/糖異生等3個核心關鍵通路，也參與了RNA降解通路，還和幾丁質酶樣蛋白一同參與了糖酵解/糖異生途徑.

對照參考組與靶標蛋白對接規(guī)則，通過PyMol軟件查看實驗組的對接結果可知，Cry1Aa，Cry4Ba，Cry1Ba，Cry5Ba，Cry1Ka，Cry1Bc，Cry1Be，Cry1Bb，Cry1Af，Cry15Aa等10個Cry蛋白可能對桃蚜具有殺蟲活性.

文中利用綜合生物信息學技術對桃蚜的關鍵抗性基因進行挖掘，并對抗蚜Cry蛋白進行預測，為后續(xù)的桃蚜關鍵抗性基因研究及新藥研發(fā)提供參考思路，具有一定的理論意義和學術意義.