捻转血矛线虫(Haemonchus contortus)是毛圆科血矛属的一种消化道线虫，主要寄生在牛、山羊、绵羊等反刍动物的皱胃。该病在世界各国广泛分布，给畜牧业造成极大的经济损失^[1-3]。近二十年来，对该病的防治主要依靠阿苯达唑等化学药物，并取得较好的效果。然而，随着该药物的频繁和不规范使用，导致捻转血矛线虫对这类药物的抗性不断增加，日益成为反刍动物消化道线虫病防治的重要问题^[4]。目前，中国、美国、加拿大、英国、法国、德国、澳大利亚、丹麦、新西兰和印度等二十多个国家均有关于捻转血矛线虫耐药性的报道，随着不同驱虫药物的使用，还存在多重耐药性流行的问题^[5]。这些报道提示多药耐药性现已成为一个全球性现象，捻转血矛线虫耐药性问题已经成为“不可忽略的真相”^[6]。因此对该病的研究应致力于探索捻转血矛线虫耐药机制的基础上，寻找有效的早期检测耐药性方法和制定新的防治措施。

据报道捻转血矛线虫对阿苯达唑产生耐药主要是由于β微管蛋白同种Ⅰ型的SNP，即F167Y、E198A、F200Y，使得阿苯达唑不能与靶点结合，从而导致耐药性的产生^[7]。然而，β微管蛋白同种Ⅰ型的SNP并不能完全解释耐药性问题，耐药性的产生可能存在多种机制。例如，Yilmaz等^[8]研究表明Cytochrome P450参与阿苯达唑耐药性的产生，Matoušková等^[9]研究表明阿苯达唑影响寄生虫酶的分泌，阻碍ATP的产生，从而导致虫体失去能量死亡，另外，还有证据表明阿苯达唑的抗药性与寄生线虫体内的Pgp基因有关，该基因可以将药物排到细胞外，从而导致耐药性的产生。迄今为止，检测阿苯达唑耐药性的方法主要包括粪便虫卵减少试验、虫卵孵化试验和体内药效试验等，但这些方法均不能在耐药性明显之前就检测出来。因此，转录组学在捻转血矛线虫耐药性研究中意义较大，本课题组前期利用Illumina HiSeq平台分析了捻转血矛线虫耐药株给药前后转录组学表达谱变化，通过GO富集与KEGG通路分析发现部分耐药相关基因^[10]。为进一步丰富耐药机制的转录组信息，本研究通过比较转录组学分析，探讨捻转血矛线虫敏感株和耐药株在代谢水平和基因功能的差异，筛选关键的耐药相关基因，旨在为深入研究捻转血矛线虫耐药机制及筛选耐药性早期检测的分子标记提供参考。

1 材料与方法 1.1 捻转血矛线虫阿苯达唑敏感和耐药虫株的采集

捻转血矛线虫阿苯达唑敏感株和耐药株均采自乌兰察布市察右后旗(北纬41°3′~41°59′，东经112°42′~113°30′)，根据粪便虫卵减少试验和对照试验采集捻转血矛线虫阿苯达唑耐药株和敏感株^[11]。

1.2 RNA提取及Illumina测序

选取敏感雌虫和雄虫，耐药雌虫和雄虫各30条，利用TRIzol方法分别对RNA进行分离和纯化，将经过质控的雌虫和雄虫核苷酸等量混合，取5 μg总RNA，使用Illumina HiSeq 4000平台进行2×150 bp测序^[12]。

1.3 测序数据分析

对测序所得的原始序列进行过滤，去除测序接头及低质量序列，得到干净序列。利用TopHat^[13]软件，将过滤后获得的干净序列比对捻转血矛线虫参考基因组，根据FPKM值对两个样本基因表达量进行归一化。采用DEGseq法^[14]，以|log2(fold change)|≥1且P≤0.05作为阈值界定敏感虫株和耐药虫株之间的差异表达基因(differentially expressed genes, DEGs)。将DEGs向Gene Ontology(GO)数据库各term进行映射，以P＜0.05为筛选条件，满足条件的term即为在DEGs中显著富集的GO term。同时将DEGs进行KEGG pathway分析，对基因涉及的显著变化的生物学调控通路进行分析。

1.4 差异表达基因的荧光定量PCR验证

选取6个差异表达基因，以β-tubulin为内参基因，采用TB Green^TMPremix Ex Taq^TMⅡ在Quant Studio 7 Flex荧光定量PCR系统中进行检测，每个样品设置3个生物学重复。

2 结果 2.1 转录组测序数据分析

对捻转血矛线虫敏感株和耐药株分别进行转录组测序，将得到的数据进行质控和组装，结果显示，敏感株和耐药株均获得8 943万以上的原始数据，且有效数据比例均在96%以上，Q30均在97%以上，表明测序质量较好。

2.2 差异表达基因筛选

以|log2(fold change)|≥1且P≤0.05为标准筛选显著性DEGs，共获得218个显著DEGs，以log2(差异倍数)为横坐标，-log10(P值)为纵坐标，对显著的DEGs绘制火山图(图 1)，结果显示，在显著DEGs中，耐药虫株中表达上调基因187个。与耐药相关的部分基因，如谷胱甘肽S转移酶、细胞色素P450、ABC转运蛋白、UDP葡萄糖转移酶、线虫角质层胶原蛋白、热休克蛋白、脂肪酶、丝氨酸苏氨酸蛋白激酶、NADH脱氢酶等9个基因的表达均在两样本中存在差异。

2.3 差异表达基因的GO功能富集分析

将测序获得的DEGs与GO数据库进行比对，注释测序基因的功能。结果显示，经GO功能富集分析后，敏感虫株和耐药虫株共有218个DEGs显著(P＜0.05)富集到206条GO terms上，其中生物学过程95条，细胞组分38条，分子功能73条，具有统计学显著意义的共63条。在BP分类中，代谢过程、黄酮类化合物葡萄糖醛酸化、黄酮生物合成过程、G蛋白偶联受体信号通路等GO terms富集水平最高；在细胞组分中，细胞外区域、细胞外间隙、细胞内膜结合细胞器等GO terms富集水平最高；在分子功能中，葡萄糖醛酸基转移酶活性、转移酶活性、转移己糖基、肌球蛋白V结合等GO terms富集水平最高(图 2)。

2.4 差异表达基因KEGG通路分析

将218个DEGs进行KEGG pathway分析，结果显示，DEGs被富集到视黄醇代谢、卟啉和叶绿素代谢、抗坏血酸和新陈代谢、细胞色素P450对异生素的代谢、药物代谢-细胞色素P450、戊糖和葡萄糖醛酸的相互转化、戊糖和葡萄糖醛酸的相互转化、药物代谢-其他酶、淀粉和蔗糖代谢、细胞凋亡-多种、溶酶体等31条通路上，其中10条显著富集的通路。

2.5 差异表达基因荧光定量PCR验证

本试验从DEGs中随机筛选6个进行荧光定量PCR验证，结果显示，经荧光定量检测的6个DEGs与RNA-Seq结果趋于一致(图 3)。

3 讨论

本试验对捻转血矛线虫阿苯达唑敏感株和耐药株进行比较转录组学分析，有效数据比例均在96%以上，Q30均在97%以上，表明测序质量良好，测序结果符合后续研究的需求。以|log2(fold change)|≥1且P≤0.05作为筛选阈值，共获得218个显著差异表达基因，并对其进行GO功能富集和KEGG信号通路分析。GO功能分析表明，差异基因主要富集在代谢过程、细胞外区域、葡萄糖醛酸基转移酶活性等，在耐药虫株中上调表达基因显著富集到代谢过程、ATP结合、GTP结合、微管细胞骨架组织、突触传播和谷氨酸调控、氯离子跨膜转运、蛋白水解、细胞凋亡、氧化还原过程等功能。其中ATP结合、GTP结合、微管细胞骨架组织、氯离子跨膜转运等功能基因都显著上调，表明在耐药虫株中耐药基因表达量上调从而导致耐药性的发生；蛋白水解、氧化还原过程和细胞凋亡等功能基因表达量也上调，可能是由于在获得耐药虫株前经过给药筛选的缘故，部分耐药基因为抵抗药物选择压力，从而导致耐药基因表达量上调。KEGG通路分析表明，显著富集的通路主要有细胞色素P450对异生素的代谢、药物代谢-细胞色素P450、药物代谢-其他酶、细胞凋亡、抗坏血酸和新陈代谢、戊糖和葡萄糖醛酸的相互转化、溶酶体等。

P450酶系是药物代谢过程中最重要的一个系统，细胞色素P450在昆虫中研究较多，主要在生长繁殖及解毒代谢过程中起着重要的作用。Zhang等^[15]对缺乏P450与过表达P450基因的寄生虫进行比较，表明后者对驱虫药表现出更强的耐药性，前人的研究结果与本试验的结果基本一致。在药物代谢-其他酶通路下包含谷胱甘肽S转移酶和UDP-糖基转移酶基因，这两个基因在生物体内主要起解毒作用，目前在螨虫等生物中均有GST基因参与耐药性的报道^[16]，本试验中该基因在耐药虫株中的表达量显著高于敏感虫株。此外，本试验中UGT基因在耐药虫株中显著上调，Vokřál等^[17]通过对捻转血矛线虫敏感和耐药虫株进行HPLC分析，结果显示，UGT基因在耐药虫株中的表达量是敏感虫株中的5倍。本试验研究结果与前人的研究结果一致，所以推测这些基因在捻转血矛线虫对阿苯达唑耐药中起到了至关重要的作用。因此，捻转血矛线虫对阿苯达唑耐药机制需进一步研究。

4 结论

通过捻转血矛线虫阿苯达唑敏感株和耐药株比较转录组学分析，获得了218显著差异表达基因，根据GO功能注释及KEGG信号通路分析，鉴定出9个耐药基因，其中ABC、GST和P450基因均在耐药虫株中差异表达；qRT-PCR验证部分差异表达基因，与转录组数据趋势一致，表明转录组数据可信度较高。研究结果为进一步探索捻转血矛线虫对阿苯达唑的耐药机制提供参考。