柞蚕微孢子虫感染后柞蚕卵转录组及免疫相关基因功能分析

徐 欣, 吴玉娇, 于 滨, 孟宪志, 陈 杰,刘中文, 张永君, 潘国庆,*

(1.西南大学, 微孢子虫感染与防控重庆市重点实验室, 重庆 400715; 2. 河南省蚕业科学研究院, 郑州 450008)

谷胱甘肽S-转移酶(glutathione S-transferases, GSTs)广泛分布于哺乳动物、植物、鸟类、昆虫、寄生虫及微生物体内,不仅能催化生物体内有害的极性化合物与谷胱甘肽(glutathione, GSH)相结合,而且也可促使体内各种具有潜在毒性的化学药物、染料及某些致癌物排出体外(Ketterer, 2001; Sheehanetal., 2001)。GSTs表达增加可能有2种机制:mRNA水平的提高和基因扩增,但基因扩增仅有少量报道(Wangetal., 1991; Zhou and Syvanen, 1997; Vontasetal., 2002)。目前已有少量确切的证据可证明GSTs的过量表达与抗性有关。如小菜蛾PlutellaxylostellaPxGST3基因所编码的酶可以降解有机磷杀虫剂,它的表达量增加与抗性有关(Huangetal., 1998)。MdGST6A在家蝇Muscadomestica的有机磷杀虫剂抗性品系中过量表达,重组酶可以与甲基对硫磷和林丹共轭(Weietal., 2001)。果蝇Drosophila的GSTd1基因在DDT抗性品系中表达量增加,重组GSTd1有DDT脱氯化氢酶活性(Tangetal., 1995)。Ranson等(2001)研究发现在冈比亚按蚊Anophelesgambiae抗性品系中Ⅲ类 GST aggst322过量表达5倍,表达的重组酶可以降解DDT。

目前哺乳动物中发现了8类细胞质GSTs 亚基(alpha, mu, pi, theta, sigma, zeta, kappa和omega)和1类微粒体GSTs亚基(马晓英等, 2007)。omega家族GST具有谷胱甘肽过氧化物酶以及巯基氧化还原酶和脱氢抗坏血酸还原酶活力,参与昆虫体内氧化胁迫修复过程,是昆虫抗氧化酶系统中的重要成员(张学尧等, 2012)。Yamamoto等(2005, 2009)对家蚕Bombyxmoriomega家族GST的分子及生化特性进行了研究,发现其参与脂质过氧化反应,并与家蚕对有机磷杀虫剂杀螟松的抗性相关,在家蚕完成外源物质解毒过程中起重要作用。

维生素C(vitamin C, VC),又称为抗坏血酸(L-ascorbic acid, AsA),在抗坏血酸氧化酶(ascorbate oxidase, AAO)和抗坏血酸过氧化酶(ascorbate peroxidase, APX)作用下脱氢被氧化成单脱氢抗坏血酸(monodehydroascorbate, MDHA)和脱氢抗坏血酸(dehydroascorbate, DHA),又分别在单脱氢抗坏血酸还原酶(monodehydroascorbate reductase, MDAR)和脱氢抗坏血酸还原酶(dehydroascorbate reductase, DHAR)作用下重新被还原成VC,使得VC得以再生(刘永立等, 2006)。其中,DHAR与MDAR一样是使抗坏血酸再生的酶,DHAR在家蚕中的预测同源基因是GST基因(杨洋等, 2016),在家蚕中的已有研究认为此基因是noppera-bo,参与家蚕蜕皮激素的合成和幼虫的发育,当此基因缺失时会导致幼虫的生长发育抑制以及角质层为光滑的表现(Enyaetal., 2015; Waltersetal., 2009)。

结合本团队感病后柞蚕Antheraeapernyi卵的转录组分析发现,柞蚕GST基因是作为感病后柞蚕体内表达差异较大的一种基因,同时其作为VC代谢通路的关键酶,在柞蚕研究中尚未被克隆。因此,本研究在家蚕GST各基因及编码的蛋白质的研究基础上,在柞蚕中克隆GST基因,并进行生物信息学分析,分析其结构和其他昆虫GSTs的进化关系,同时检测柞蚕GST基因在柞蚕微孢子虫侵染后柞蚕产卵0-10 d时卵中的表达量变化,获得柞蚕GST基因的基本特性。

1 材料与方法

1.1 供试材料

柞蚕品种为豫大一号,由河南省蚕业科学研究院繁育提供。

1.2 主要试剂和仪器

动物组织总RNA提取试剂盒RNAprep Pure Tissue Kit,FastQuant cDNA第1链合成试剂盒FastQuant RT Kit(with gDNase),质粒提取试剂盒,凝胶回收试剂盒,2×Taq PCR MasterMix,FastFire 快速荧光定量 PCR 预混试剂(SYBR Green),琼脂糖,GeneGreen核酸染料均为天根生化科技(北京)有限公司产品。试验中所需的限制性核酸内切酶、连接酶、高保真酶、大肠杆菌Escherichiacoli感受态细胞DH5α、载体 pMD-18T、pET28a等均购自TaKaRa宝日医生物技术(北京)有限公司。引物合成与基因测序工作委托生工生物工程(上海)股份有限公司完成。大肠杆菌用培养基及其他常用试剂为国产或进口分析纯。

TGL16M高速冷冻离心机,购自湖南湘立科学仪器有限公司;PCR仪(9200型),购自Applied Biosystems公司;超微量分光光度计,购自美国Thermo Scientific公司;1600凝胶成像仪,购自上海天能科技有限公司;超声波破碎仪,购自宁波新芝生物科技股份有限公司;电泳电源、电泳槽、水平摇床,购自北京市六一仪器厂。

1.3 感染柞蚕微孢子虫后柞蚕卵转录组及数据分析

取同一批柞蚕蛹,随机选取镜检携带柞蚕微孢子虫的母蛾10头,未感染柞蚕微孢子虫正常母蛾10头作为对照,待产卵后,将柞蚕卵进行表面消毒,从中各随机选取柞蚕卵200粒作为1个重复,每组3个重复,放入液氮中保存备用。用DEPC水把研钵、剪刀、镊子等用具洗干净,放在恒温干燥箱内,160 ℃烘烤4 h备用。再使用动物组织总RNA提取试剂盒,按照其说明书方法从正常对照柞蚕卵和感染柞蚕卵中分别提取总RNA,置于-80 ℃中保存备用。

用mRNA富集法或rRNA去除法对总RNA进行处理(mRNA富集:用带有OligodT的磁珠富集有polyA尾巴的mRNA;rRNA去除:用DNA探针杂交rRNA),RNaseH选择性消化DNA/RNA杂交链,再用DNaseⅠ消化掉DNA探针,纯化后即得到所需RNA;用打断Buffer把获得的RNA片段化,随机的N6引物进行反转录,再合成cDNA二链形成双链DNA;把合成的双链DNA末端补平并5′端磷酸化,3′端形成突出一个“A”的粘末端,再连接一个3′端有凸出“T”的鼓泡状的接头;连接产物通过特异的引物进行PCR扩增;PCR产物热变性成单链,再用一段桥式引物将单链DNA环化得到单链环状DNA文库;上机测序。

过滤掉低质量、接头污染以及未知碱基N含量过高的原始读段,过滤后获得有效读段。然后将有效读段比对到柞蚕参考基因组上进行新转录本预测、SNP &InDel和差异剪接基因检测。得到新转录本之后,将具有蛋白编码潜力的新转录本加入到柞蚕参考基因组序列中构成一个完整的参考序列,然后计算基因表达量。最后,对于多个样品根据需求检测不同样品之间的差异表达基因,并对差异表达基因做深入的聚类分析和功能富集分析。

根据GO功能注释结果以及官方分类,将差异表达基因进行功能分类,同时使用R软件中的phyper函数进行富集分析,计算P值,然后对P值进行 FDR校正, 校正后的P值又称Q值,通常Q≤0.05的功能视为显著富集,以Q值最小的前20个GO功能条目作图。富集气泡图从3个维度展示GO功能条目的富集程度,默认以Q值最小的前20个GO功能条目或按照选择的GO功能条目(按Q值排序,最多60个)作图。

根据KEGG注释结果以及官方分类,将差异表达基因进行生物通路分类,其余处理同GO功能注释。

1.4 免疫相关基因ApGSTo1的克隆及生物信息学分析

根据1.3节转录组中差异表达基因分析结果,选择保守结构区域,结合PCR引物设计原则,用Primer Premier 5.0设计ApGSTo1特异引物(表1),并选择ApA(AntheraeapernyiActin)为内参基因(表1),其退火温度为50 ℃,循环数为30。

表1 引物信息Table 1 Primer information

使用动物组织总RNA提取试剂盒,按照其说明书方法从感染柞蚕微孢子虫的柞蚕蛹中提取总RNA,按照FastQuant cDNA第一链合成试剂盒说明书的步骤将RNA反转为cDNA, PCR反应体系: 2×Taq PCR Mastermix 12.5 μL, 正反向引物(10 μmol/L) 各1 μL, cDNA 1 μL, ddH2O加至25 μL。反应程序: 94 ℃预变性10 min; 94 ℃变性30 s, 50 ℃退火30 s, 72 ℃延伸30 s, 共35个循环;72 ℃再延伸10 min。反应结束后用1.0%琼脂糖凝胶电泳检测PCR产物,回收目的片段,回收产物与pMD18-T载体在16 ℃过夜连接,连接产物转化大肠杆菌DH5α,在含氨苄青霉素的培养基上进行蓝白斑筛选,随机挑取白色菌落进行菌落鉴定,筛选出阳性克隆,由上海生工生物工程技术服务有限公司进行测序鉴定。

柞蚕ApGSTo1信号肽及跨膜区分析使用在线软件SMART(http:∥smart.embl-heidelberg.de/),并运用PSIPRED网站(http:∥bioinf.cs.ucl.ac.uk/psipre.html)进行二级结构在线预测。蛋白质高级结构模型预测利用在线软件SWISS-MODEL(https: ∥swissmodel.expasy.org/)。为研究柞蚕ApGSTo1蛋白与在鳞翅目不同昆虫omega家族GST之间的进化关系,本研究利用果蝇基因组数据库(FlyBase)、家蚕基因组数据库(Silkworm Genome Database),运用SnapGene软件和NCBI等数据库,选取已报道的鳞翅目不同昆虫GSTo1的氨基酸序列,与柞蚕ApGSTo1氨基酸序列用SnapGene软件进行多序列比对,并用MEGA 11.0软件分析系统进化关系。

1.5 感染柞蚕微孢子虫柞蚕卵不同发育阶段ApGSTo1表达量的RT-qPCR检测

随机选取镜检感染柞蚕微孢子虫母蛾20头,未感染的正常母蛾20头,待产卵后分成感染组和正常对照组两组,置于室温25 ℃的柞蚕卵孵化室进行正常孵化发育,在柞蚕卵从母蛾体内产出的0-10 d的时间内,每天从两个组中随机选取蚕卵50粒,将蚕卵进行表面消毒,放入液氮中保存备用。待10 d两个组的样品取完后,提取总RNA,并反转录为cDNA,方法同1.3节。

根据克隆测序得到的ApGSTo1序列,选择保守结构区域,结合RT-qPCR引物设计原则,用Primer Premier 5.0设计ApGSTo1的RT-qPCR特异引物(表1),目的片段长度为223 bp;ApA为内参基因,设计其qRT-PCR特异引物(表1),目的片段长度为134 bp。根据FastFire 快速荧光定量 PCR 预混试剂(SYBR Green)的实验方法进行qRT-PCR。反应程序: 95 ℃ 1 min; 95 ℃ 5 s, 60 ℃ 10 s, 72 ℃ 15 s, 循环数为40。

1.6 数据分析

采用 2-ΔΔCt法进行基因相对表达量分析。采用 GraphPad Prism 软件进行数据处理和作图,并用SPSS 26.0数据分析软件单因素方差分析(one-way ANOVA)进行显著性分析。

2 结果

2.1 感染柞蚕微孢子虫后柞蚕卵转录组差异表达基因的GO功能注释

将感染柞蚕微孢子虫与未感染的正常对照组的柞蚕卵进行对比,每组3个重复,在6个样品的转录组测序结果中,预测的新基因为3 400个;共检测到的表达基因数为17 453,其中已知的基因为 14 152个,预测的新基因为3 301个。其中,在转录本的长度分布分析中发现超过3 000 nt序列长度的基因数量最多,将近3 500个,其他长度区间的基因数量在1 000以下,分布较均匀。同时,根据各个样品基因表达水平,检测到的显著差异表达基因数目为89个。对每个比较组差异基因的FPKM值进行聚类,用热图显示如图1。

图1 感染柞蚕微孢子虫后柞蚕卵转录组中差异表达基因表达聚类热图Fig. 1 Expression clustering heatmap of differentiallyexpressed genes in the transcriptome of Antheraea pernyieggs after infection with Nosema pernyiLW1A, LW2A, LW3B: 感染柞蚕微孢子虫后柞蚕卵组The groups of A. pernyi eggs infected with N. pernyi.

将感染柞蚕微孢子虫后的柞蚕卵与正常对照组的柞蚕卵差异表达基因进行GO功能注释,在生物学过程(biological process)类别中注释到细胞过程(cellular process)和代谢过程(metabolic process)的基因数量最多;细胞组分(cellular component) 类别中注释到细胞解剖实体(cellular anatomical entity)的基因数量最多;分子功能(molecular function) 类别中注释到催化活性(catalytic activity)的基因数量最多,而且是所有分类中基因数量最多的一项,约19个(图2)。

图2 感染柞蚕微孢子虫后柞蚕卵转录组中差异表达基因的GO功能注释Fig. 2 GO functional annotation of differentially expressed genes in the transcriptome of Antheraea pernyi eggs after infection with Nosema pernyi

从差异表达基因GO富集气泡图(图3)中看出,显著性较高的差异表达基因大多数属于各种代谢过程(metabolic process),但是富集因子数值小,说明在通路中富集显著性的可靠度较低。富集到细胞外基质结构成分(extracellular matrix structural constituent)和胶原三聚体(collagen trimer)的基因数量也较多,显著性和富集因子数值也较高。

图3 感染柞蚕微孢子虫后柞蚕卵转录组中差异表达基因GO富集气泡Fig. 3 GO enrichment bubble diagram of differentially expressed genes in the transcriptome of Antheraea pernyi eggs after infection with Nosema pernyi

2.2 感染柞蚕微孢子虫后柞蚕卵转录组差异表达基因的KEGG注释

根据KEGG注释结果以及官方分类,将感染柞蚕微孢子虫后的柞蚕卵与正常对照组的柞蚕卵转录组中差异表达的基因进行生物通路分类,分为细胞过程(cellular processes)、环境信息处理(environmental information processing)、遗传信息处理(genetic information processing)、代谢(metabolism)和生物体系统(organismal systems)共5个分支。细胞过程的通路中参与运输和分解代谢(transport and catabolism)的基因数量最多,环境信息处理的通路中参与信号转导(signal transduction)的基因最多,遗传信息处理的通路中参与翻译(translation)的基因最多,代谢中参与全局和概述图谱(global and overview maps)的基因数量最多,生物体系统中参与消化系统(digestive system)的基因最多。与GO侧重于单个基因的功能分析相比,KEGG分析则侧重于寻找和分析生物体系统中基因之间的相互作用,在感染柞蚕微孢子虫后差异基因的KEGG分析中,参与运输和分解代谢、信号转导和消化系统这3个通路的基因数量都达到12个,都是后期重点筛选的通路。

KEGG通路的富集气泡图(图5)中,生热作用(thermogenesis)、蛋白质消化和吸收(protein digestion and absorption)两个通路聚集的基因数目较多,显著性也高。结合差异基因表达量的结果,log2FC值为7.501,根据KEGG的id号,匹配到的为Tudor结构域蛋白;根据表达量差异值的排序,依次匹配的为:MFS转运蛋白、RB1诱导型卷曲螺旋蛋白、溶质载体PCFT/HCP系列蛋白、醛脱氢酶、溶酶体酸性脂肪酶/胆固醇酯水解酶、WD重复和FYVE结构域蛋白、BRG1相关因子、谷胱甘肽S-转移酶,其中GSTs排第12位,log2FC值为3.979,主要参与谷胱甘肽代谢过程。然后,在差异最显著的基因中寻找与昆虫抗逆相关的目标基因,发现GST基因在昆虫和家蚕中均有参与氧化胁迫修复过程、抗性相关、完成外源物质解毒作用过程等相关研究结果,因此将GST基因作为目标基因进行下一步试验。

图5 感染柞蚕微孢子虫后柞蚕卵转录组中差异表达基因 KEGG通路富集气泡Fig. 5 Enrichment bubble diagram of KEGG pathways of differentially expressed genes in the transcriptome of Antheraea pernyi eggs infected with Nosema pernyi

2.3 ApGSTo1的克隆与序列

PCR获得ApGSTo1的cDNA序列(GenBank登录号: OP661169.1),编码区全长767 bp,编码254个氨基酸。在SMART程序中搜索到ApGSTo1的5个蛋白结构域(GST_N, GST_N_3, GST_N_2, GST_C_3和GST_C_2),这5种组合形成6种柞蚕ApGSTo1的预测结构,其中2种蛋白结构域组成为:GST_N(第17-90位氨基酸)和GST_C_2(第121-205位氨基酸)。还有GST_N_3, GST_N_2, GST_C_3这3种蛋白结构域。同时还有未显示的其他6种蛋白结构域(Glutaredoxin, Tom37, GST_N_4, GST_C, GST_C_6和Tom37),以及3种SMART区域(Amelogenin, GIT和MHC_II_alpha)。

在PSIPRED网站上柞蚕GSTo1的二级结构预测结果显示,该序列中存在12个α螺旋和4个β折叠,在 N 端存在1个β-α-β-α-β-β-α的结构基元,如图6所示。

图6 柞蚕ApGSTo1的氨基酸序列及二级结构预测Fig. 6 Prediction of amino acid sequence and secondary structure of ApGSTo1 of Antheraea pernyi

将柞蚕ApGSTo1序列提交到SWISS-MODEL服务器构建模拟结构,系统对该模型的QMEANscore4评分值为为0.78,如图7(A)。模拟结构中柞蚕GSTo1与家蚕omega家族GST晶体结构(图7: A, B)相似性为41.95%, 与人omega家族GST晶体结构(图7: A, C)相似性为33.77%,结果显示其拥有典型的omega家族GST特征,N端结构域由4个β折叠和3个α螺旋组成,C端结构域由8个α螺旋组成,二级结构元件相互之间以转角和无规则卷曲相互连接,这表明柞蚕GSTo1在蛋白整体结构上与omega家族GST结构类似。

图7 柞蚕ApGSTo1的蛋白质结构预测Fig. 7 Protein structure prediction of ApGSTo1 of Antheraea pernyiA: 柞蚕GSTo1晶体结构GSTo1 crystal structure of A. pernyi; B: 家蚕omega家族的GST晶体结构Omega-class GST crystal structure of Bombyx mori; C: 人omega 家族的GST晶体结构Omega-class GST crystal structure of Homo sapiens.

由图8可知,所有序列聚为3个大类,柞蚕ApGSTo1(GenBank登录号: OP661169.1)聚类在GSTo1一组,证明ApGSTo1属于omega 1家族GST,并与家蚕GSTo1被聚到一个分支中,表明亲缘关系最近;在GSTo1分组中,从遗传距离来看,鳞翅目夜蛾科的甜菜夜蛾Spodopteraexigua、斜纹夜蛾S.litura和草地贪夜蛾S.frugiperda的 GSTo1遗传聚类最远;其他两个Omega家族的序列分别聚类为GSTo2和GSTo3组,且GSTo2与GSTo1组分为一个支,GSTo3单独分支,说明GSTo1与GSTo2家族的亲缘关系较近。

图8 基于氨基酸序列构建的柞蚕ApGSTo1与不同鳞翅目昆虫omega家族的GSTs的系统进化树Fig. 8 Phylogenetic tree of Antheraea pernyi ApGSTo1 and omega-class GSTs of different lepidopteran insects based on amino acid sequence

2.4 ApGSTo1在柞蚕卵不同发育时期的表达变化

在柞蚕卵从母蛾体内产出的0-10 d内,ApGSTo1的表达量在正常对照组和感染柞蚕微孢子虫组内的整体变化趋势都是先升高再降低(图9),都是在产卵后3 d时达到最高值,这应该与胚胎的发育规律有关;在产卵后2-6 d时,感染柞蚕微孢子虫组的基因表达量比正常对照组的都高,而且在产卵后2 d时与正常对照组差异极显著(P<0.01),在产卵后6 d时显著高于正常对照组的(P<0.05),这应该与ApGSTo1作为巯基转移酶和抗坏血酸氧化还原酶参与柞蚕卵内的抗逆胁迫反应有关,产卵后2 d时也是柞蚕卵发育与抗逆的关键时期。

3 讨论

柞蚕微粒子病是生产上的重要病害,威胁种质资源的保育及生产的安全,因此对于柞蚕微孢子虫分离检测、宿主应答、传播途径、综合防治等方面的研究一直是热门方向,转录组学和蛋白质组学的研究,为探究柞蚕微孢子虫的免疫应答模式、柞蚕微粒子病的发病机理和开发防治药剂提供直接依据。本研究首次获得了柞蚕微孢子虫感染柞蚕卵与未感染柞蚕卵的转录组数据,检测到的表达基因数为17 453,其中已知的基因为 14 152个,预测的新基因为3 301个,补充了柞蚕在分子水平上的抗病性研究的空白,再结合柞蚕放养习性和家蚕的室内饲养的差别,柞蚕在抗病性研究领域将是更理想的模式昆虫,因此后续除了对已知基因的挖掘外,还将对新基因进行结构和功能的预测,以期发现更多的抗病性相关通路或蛋白。

昆虫的解毒酶系主要包括: 谷胱甘肽-S-转移酶系、细胞色素P450 酶系和酯酶(如羧酸酯酶、乙酰胆碱酯酶)等(Kerkut and Gilbert, 1985)。其中,谷胱甘肽-S-转移酶是广泛分布于动物、植物、酵母和细菌等生物体中的一类由多基因编码的、具有多种功能的重要酶系。在转录组GO功能注释中发现分子功能类别中的差异表达基因数量最多,其中参与催化活性的基因数量达19个(图2),再结合KEGG分析中差异基因数量最多的3个通路(运输和分解代谢、信号转导和消化系统)的结果(图4)可知,感染柞蚕微孢子虫后的柞蚕卵发生的主要变化集中在物质代谢上,与免疫或者抗逆性直接相关的应该是重要代谢通路的酶的表达量变化。同时根据差异表达基因log2FC值大小顺序,差异最显著的一些基因参与的KEGG通路也多数为转运蛋白、代谢酶、结构蛋白等,如MFS转运蛋白、醛脱氢酶、溶酶体酸性脂肪酶和谷胱甘肽S-转移酶;再结合前人关于昆虫的解毒酶系的研究结果发现,最可能参与柞蚕免疫相关的候选基因为谷胱甘肽S-转移酶基因,而且有研究证明了在昆虫中,GSTs的主要功能是催化谷胱甘肽的巯基与一些亲电子类有毒物质进行轭合反应,从而保护一些蛋白质免受损伤,达到解毒目的(赵国栋等, 2010; 张婷等, 2011)。

家蚕微孢子虫发育分为感染期、裂殖增殖期和孢子形成期3个阶段(Ishiharaetal., 1969)。前人猜测卵龄20 h的带毒蚕卵内病原体主要以裂殖体为主,王璐等(2019)采用实时荧光定量 PCR分别对阳性对照、常规即时浸酸和高温即时浸酸处理后卵龄 24-240 h蚕卵中的病原数量逐日检测,发现阳性对照组卵龄48 h蚕卵中的病原量出现急剧增长,推测产卵后 24-48 h 是家蚕微孢子虫的增殖高峰期,此后每日病原数量无显著变化,可能与蚕卵逐渐进入滞育状态,病原生长发育同样受到抑制有关。这与本研究中ApGSTo1的表达量在正常对照组和感染柞蚕微孢子虫组内的整体变化趋势(图9)一致,同时发现在柞蚕卵从母蛾体内产出的0-10 d内,基因表达量在2 d时极显著高于正常对照组,在6 d时显著高于正常对照组(图9),推测产卵后2 d也是柞蚕卵发育与抗逆的关键时期,这也与家蚕微孢子虫增殖高峰期的时间一致。寻找和培育抗微粒子病新品种一直是柞蚕育种上的一大难关,由于缺乏靶标检测基因,无法实现从分子育种角度对现有柞蚕品种进行筛选,因此ApGSTo1可以作为柞蚕幼虫免疫柞蚕微孢子虫的靶标基因,为后续寻找和培育抗微粒子病新品种提供帮助。