桃蚜关键抗性基因挖掘及抗蚜Cry蛋白预测

瞿贵军，林毅

(华侨大学 化工学院，福建 厦门 361021)

桃蚜(Myzuspersicae)的寄主植物多达50多科400余种，桃蚜主要以果树、蔬菜、烟草、花卉等数百种植物为食，是世界性的重要经济害虫之一[1]，也是世界传播病毒昆虫之首[2].桃蚜的次级共生菌主要是保护其宿主免受病原体和天敌的侵害，增强对杀虫剂的抗性、新陈代谢和生物合成等[3].桃蚜的生命周期复杂，目前主要通过化学杀虫剂杀桃蚜，但桃蚜对使用过的杀虫剂都能产生不同程度的耐药性，具有很强的遗传变异性和适应性[4].目前，桃蚜在抗性克隆中通过细胞色素P450(P450)和UDP葡萄糖醛酸转移酶(UGT)高度上调表达[5].环境友好型生物农药Cry蛋白已被国际抗性行动委员会(IRAC)分类方案收录，作为杀虫剂活性成分[6]，在生物农药研发中，Cry41-related蛋白可通过增强组织蛋白酶B活性杀死桃蚜[7]，而Cry1Cb2蛋白也有望成为一种有前途的抗蚜蛋白[8].

近年来，在基因挖掘和调控机制等方面的研究中，加权基因共表达网络分析(WGCNA)[9]、差异表达基因分析(DEGs)[10]、GO/KEGG富集分析、蛋白质互作网络分析(PPI)等综合生物信息技术分析取得了许多研究成果.He等[11]通过WGCNA揭示昆虫真菌病原体的各种不利胁迫适应性响应的基因和信号通路.Li等[12]通过差异表达基因分析，为评估传毒和非传毒桃蚜对植物病毒的反应机制提供依据.Chin等[13]使用Cytoscape软件中的CytoHubba插件，根据网络特征对网络中的节点进行排名.基于此，本文对GSE18659，GSE37310等2个GEO数据(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数据(16个样本)进行在线分析，离群样本阈值设定为0.00，共表达网络类型选择unsign，hub相关性阈值设定为0.9，网络权重阈值设定为0.00，模块最小基因数设定为20.00，模块合并阈值设定为0.25，敏感性(1～4)设定为2.00，软阀值设定为0.00，获得关键抗性基因文件和基因互作文件.

然后，采用Cytoscape软件的Cytohubba插件，对已用R语言筛选出weight值大于0.48的基因互作文件进行打分，通过MCC算法和Degree算法，得到两种算法各前10名基因探针对应的基因.

最后，通过GEO数据库自带基于limma包R语言的GEO2R工具，对GSE18659进行差异表达基因分析，根据校正后的P值Padj<0.05，差异倍数|log FC| >1，筛选出差异表达的关键抗性基因.在DAVID网站(https:∥david.ncifcrf.gov)选择桃蚜研究的模式生物豌豆蚜Buchnera aphidicola str.APS (Acyrthosiphonpisum)生物体进行在线分析，获得GO/KEGG富集分析数据.

1.3 抗性调控网络的构建

将挖掘到的关键抗性基因通过String数据库(https:∥string-db.org)构建抗性调控网络.将模式生物豌豆蚜(Acyrthosiphonpisum)的蛋白质互作数据(String数据库)与桃蚜数据(Uniprot数据库(https:∥www.uniprot.org))进行合并，在豌豆蚜中去掉与桃蚜名称相同的蛋白，再与WGCNA和DEGs获得的桃蚜关键抗性基因进行名称比对，根据与关键基因描述名称是一致的原则进行筛选，获得蛋白质互作网络分析所需的蛋白名.

在String数据库中进行蛋白质互作分析时，分别选择桃蚜共生菌Buchnera aphidicola str.USDA(Myzuspersicae)生物体和豌豆蚜(Acyrthosiphonpisum)生物体，通过桃蚜和豌豆蚜的相同蛋白名互换基因名实现两个调控网络串联，合并绘制桃蚜抗性调控网络.同时，通过NCBI网站提供的blastx工具，将前期WGCNA和DEGs得到的关键抗性基因序列转化为蛋白序列，在String数据库选择桃蚜共生菌Buchnera aphidicola str.USDA(Myzuspersicae)生物体进行蛋白质互作分析.将经Multiple proteins和Multiple sequences检索策略得到的蛋白互作数据合并，得到string_interactions.tsv文件.最后，采用Cytoscape软件打开该数据重新绘制，获得桃蚜抗性调控网络图；采用基迪奥生物云工具(https:∥www.omicshare.com)绘制抗性调控权重网络图，连通性方法选择硬阀值，权重范围选择默认值(0.5～1.0).

1.4 关键靶标蛋白库的构建

首先，在桃蚜抗性调控网络中筛选出DEGs中上调表达的关键基因.

然后，通过Cytoscape软件的ClueGO，CluePedia插件，选择模式生物豌豆蚜生物体对桃蚜关键抗性基因进行分析，得到值得关注的蛋白.

最后，整合蛋白互作数据的分析结论，绘制Cry蛋白杀桃蚜的靶标蛋白调控网络.采用NCBI网站的Blastp工具，在UniProt网站输入靶标蛋白名或登录号，得到同源序列.利用Swiss-MODEL网站(https:∥swissmodel.expasy.org)进行同源建模，选择“Experimental Structures”的Range覆盖率最高或“Homology model”的QMEANDisCo分数最高，确定靶标蛋白库的分子对接实验模型.

1.5 Cry抗蚜蛋白的预测

在Molprobity网站(http:∥molprobity.biochem.duke.edu/index.php)对同源建模的3D模型结果进行评估，在Z-dock网站(https:∥zdock.umassmed.edu)将评估合格的靶标蛋白与Cry蛋白的3D结构进行在线分子对接，所有参数设置为默认值，将实验结果文件导入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在某些情况下可能存在协同表达，UGT与GST则存在协同表达；在cyan模块中发现了Cathepsin B；在blue模块中发现了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模块中发现了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个上调表达基因中发现了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富集分析的数据.与桃蚜共生菌相关的关键抗性基因GO富集分析，如表3所示.表3中：n为基因数；η为基因占比.由表3可知：成功转化了58个关键抗性基因，在生物学过程中，关键抗性基因参与细胞氧化还原稳态和糖酵解过程；在细胞定位中，关键抗性基因富集于细胞质；从分子功能上看，关键抗性基因具有使ATP结合、氨基酰基tRNA编辑活性和核酸结合的作用.由于这些关键抗性基因与桃蚜抗性相关，从基因数和基因占比来看，其主要富集于细胞质和ATP结合.

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

与桃蚜共生菌相关的关键抗性基因KEGG富集分析，如表4所示.由表4可知：这些关键抗性基因主要富集于氨酰tRNA生物合成、碳代谢、不同环境中的微生物代谢、糖酵解/糖异生、RNA降解、谷胱甘肽代谢、磷酸戊糖途径、柠檬酸循环(TCA循环)等通路；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蛋白能与其发生连接.然后，将这10个Cry蛋白与Actin，APN，PFKA，Cathepsin B，Dead-CshA，DNAK等7个靶标蛋白在Z-dock网站进行分子对接实验，参考抗蚜Cry1Cb2蛋白与这7个靶标蛋白的参考组对接规则，用PyMol软件查看实验组的对接结果，均能发生连接.

3 讨论与分析

桃蚜抗性机制可以分为靶标敏感性降低和抗性代谢能力增强等两大类，通过其羧酸酯酶的过量表达、乙酰胆碱酯酶的突变、电压门控钠通道突变、GABA受体亚基基因的复制、细胞色素P450CYP6CY3的过表达、烟碱乙酰胆碱受体(NAChR)的突变和ABC转运蛋白等发挥蚜虫的抗性机制[16-22]，文中对已报道的桃蚜关键抗性蛋白UGT 和GST预测存在协同表达，P450和CarE在某些通路下协同表达，DDX52，Cathepsin B，PFKA，Hsp60 protein等与UGT 和GST可能存在协同表达，数据显示P450，CarE，GST，HSP60，DDX52，Cathepsin B等均以上调表达参与抗性.GO富集分析表明，关键抗性基因富集在细胞氧化还原稳态和糖酵解过程，从基因数和基因占比来看，其主要富集于细胞质和ATP结合，说明桃蚜共生菌参与抗性主要是在细胞质，通过与ATP结合相关的基因调控、代谢通路或信号转导等来实现.KEGG富集分析表明，eno，pfka等基因参与RNA降解通路，而eno和acpp直接与pfka互作，dead能直接作用于dnak.

除已报道的UGT 和P450等桃蚜关键抗性蛋白，在基于String数据库预测构建的桃蚜抗性调控网络中，烯醇化酶(Enolase)，PFKA，CLP，DNAK，CshA，Pescadillo，Autophagy protein和ACP等蛋白也值得重点关注.烯醇化酶是由蚜虫ERV畸胎细胞释放，介导参与病原体和肿瘤细胞侵入组织及逃避宿主免疫反应的酶的激活[23].Ae-ENO是由胚胎膜分离的细胞产生的，这些细胞和毒液一起在宿主蚜虫的调节和营养开发中发挥重要作用.PFKA是糖酵解最重要的调节酶，它和LOC100162014一同参与了糖酵解/糖异生途径，LOC100162014是模式生物豌豆蚜几丁质酶样蛋白(CLP)的基因[24].此外，Cry-related蛋白能够与桃蚜共生菌PFKA结合.DNAK是豌豆蚜细胞内共生热休克蛋白同系物[25]，DNAK伴侣蛋白也是蚜虫蜜露的蛋白质之一[26]，它可能在植物-蚜虫相互作用中充当介质.CshA是一种二聚体 DEAD-box解旋酶，可与核糖核酸酶合作进行mRNA 转换[27].上调表达的acypi006185-pa是模式生物豌豆蚜的Pescadillo homolog基因，Pescadillo是一种在恶性细胞中异常表达的细胞周期调节蛋白，在细胞增殖中起着至关重要的作用，并且可能是致癌转化和肿瘤进展所必需的[28]，研究发现DEAD-box解旋酶DDX27可以与Pes1特异性相互作用[29].模式生物豌豆蚜基因10016879在String数据库注释是自噬蛋白，参与自噬囊泡的形成，比较转录组分析高粱甘蔗桃蚜抗性的遗传机制发现在SCA感染后抗性基因型中自噬的基因上调表达[30]，另有研究发现，选择性自噬通过NBR1介导的病毒衣壳蛋白和颗粒靶向抑制花椰菜花叶病毒侵染，病毒触发的自噬很大程度上独立于 NBR1 的方式，防止受感染植物的广泛衰老和组织死亡，这种生存功能显著延长了病毒生存的时间，从而增加了蚜虫载体和 CaMV 传播获得病毒颗粒的机会[31].Autophagy protein极有可能参与了桃蚜的抗性作用.acpp是酰基载体蛋白ACP的基因，ACP是一类具有保守的丝氨酸残基的小分子量的酸性蛋白质，它也可作为抗菌剂的靶点[32]，参与桃蚜的代谢途径和次级代谢产物的生物合成途径，也有研究发现delta914:0-酰基载体蛋白脂肪酸去饱和酶基因的表达对于在抗虫天竺葵(Pelargoniumxhortorum)中发现的omega 5漆树酸是必需的[33]，参考受到广泛认可的Bravo穿孔模型[34]，Cry蛋白由于高pH值和还原条件而溶解于中肠腔，并被肠道蛋白酶激活，产生毒素片段，但目前多数Cry蛋白杀桃蚜效果并不理想，这也可能与上调表达的桃蚜acpp基因有关系.ACP和PFKA存在蛋白质互作，而PFKA参与了碳代谢、微生物代谢和糖酵解/糖异生等3个核心关键通路，也参与了RNA降解通路，还和几丁质酶样蛋白一同参与了糖酵解/糖异生途径.

对照参考组与靶标蛋白对接规则，通过PyMol软件查看实验组的对接结果可知，Cry1Aa，Cry4Ba，Cry1Ba，Cry5Ba，Cry1Ka，Cry1Bc，Cry1Be，Cry1Bb，Cry1Af，Cry15Aa等10个Cry蛋白可能对桃蚜具有杀虫活性.

文中利用综合生物信息学技术对桃蚜的关键抗性基因进行挖掘，并对抗蚜Cry蛋白进行预测，为后续的桃蚜关键抗性基因研究及新药研发提供参考思路，具有一定的理论意义和学术意义.