2. 内蒙古自治区综合疾病预防控制中心, 呼和浩特 010018;
3. 内蒙古农业大学生命科学院, 呼和浩特 010018
2. Inner Mongolia Center for Disease Control and Prevention, Hohhot 010018, China;
3. College of Life Sciences, Inner Mongolia Agricultural University, Hohhot 010018, China
捻转血矛线虫(Haemonchus contortus)是一种重要的胃肠道线虫，寄生在牛、羊等反刍动物的皱胃，可引起宿主贫血、腹泻，严重的可导致动物，尤其是幼畜大批量死亡，给全世界畜牧业带来了巨大的经济损失[1]。目前，对该病的防治主要依靠阿苯达唑、伊维菌素等化学药物，阿苯达唑于1979年在中国研制成功，由于其具有高效、广谱、低毒的特点，被广泛用于人畜寄生虫病的治疗。但几十年以来阿苯达唑的反复频繁使用和盲目用药，导致了虫体逐渐产生耐药性。因此，明确捻转血矛线虫对阿苯达唑的抗性机制对合理控制捻转血矛线虫及其抗药性发展具有重要意义。
1986年，南非首次报道寄生虫对阿苯达唑的耐药性，随后中国、美国、加拿大、英国、法国、德国、澳大利亚、丹麦、新西兰、印度等二十余个国家均相继报道[2]。近年来，对捻转血矛线虫阿苯达唑耐药性的研究多集中在粪便虫卵减少试验(faecal egg reduction，FECRT)、虫卵孵化试验(egg hatch assays，EHAs)和幼虫发育试验等传统的方法。Barrere等[3]以粪便虫卵减少试验和幼虫发育试验评估加拿大地区16个牧场绵羊胃肠道线虫抗性，结果表明，所有牧场的绵羊捻转血矛线虫均对阿苯达唑产生耐药性。Mickiewicz等[4]以粪便卵数减少试验和虫卵孵化试验对波兰山羊胃肠道线虫进行评估，结果表明，捻转血矛线虫对阿苯达唑产生了耐药性，这也是波兰首次报道山羊对阿苯达唑驱虫药的耐药性。此外，耐药性研究发现，捻转血矛线虫对阿苯达唑耐药性与β-微管蛋白同种型Ⅰ基因的单核苷酸多态性(single nucleotide polymorphisms，SNPs)相关，即E198A(GAA至GCA)、F167Y(TTC至TAC)和F200Y(TTC至TAC)的SNP，三个氨基酸位点的突变可用作抗药性检测的标记[5]。Zhang等[6]建立了ARMS-PCR检测捻转血矛线虫β-微管蛋白同种型Ⅰ基因的SNP-E198A，从而鉴别捻转血矛线虫对阿苯达唑的抗性。薄新文和李祥瑞[7]基于β-微管蛋白同种型Ⅰ基因SNPs建立了能检测捻转血矛线虫阿苯达唑药物抗性的多重PCR方法，并用虫卵孵化法从分子生物学角度证实其可行性。但这些方法均不能在耐药性明显出现前检测出来。
前人在捻转血矛线虫耐药性研究中取得了阶段性成果，但有关捻转血矛线虫的组学研究还很匮乏[8-9]，尤其是对其耐药性相关的转录组学研究仍属未知。本研究以捻转血矛线虫阿苯达唑耐药株为试验材料，比较捻转血矛线虫给药前(before albendazole administration，bBZ)和给药后(after albendazole administration，aBZ)的转录组学差异，从转录组水平研究捻转血矛线虫对阿苯达唑的抗性机制，筛选可对耐药性进行早期检测的分子标记，进而为建立线虫耐药性早期快速检测方法奠定基础。
1 材料与方法 1.1 捻转血矛线虫耐药虫株的采集捻转血矛线虫耐药虫株采自内蒙古鄂尔多斯市乌审旗绵羊皱胃中，根据雌虫和雄虫的形态特征，在体视镜下鉴定，标记好给药前、给药后及雌雄虫存入液氮中，用于RNA的提取。
1.2 RNA提取及检测选取捻转血矛线虫耐药虫株给药前雌虫、给药前雄虫、给药后雌虫及给药后雄虫各1管，每管60条，参照Qiagen+TRK1002(LC Science，Houston，TX)说明书对给药前后的捻转血矛线虫雌、雄虫分开提取总RNA，使用RQ1酶消解RNA中的DNA，用Bioanalyzer 2100和RNA 6000 Nano Lab Chip Kit(Agilent，CA，USA)进行总RNA质量和浓度检测。
1.3 去rRNA链特异性文库构建经过质检合格的总RNA(RIN>7.1)，用Epicentre Ribo-Zero Gold试剂盒去除核糖体RNA(rRNA)，纯化回收剩余部分rRNA depleted RNA，并将其随机打断成短片段，以此为模板，合成cDNA。对反转录获得的cDNA进行PCR扩增后，采用Illumina Hiseq4000测序平台进行文库的测序，测序由杭州联川生物技术股份有限公司完成。
1.4 测序数据分析首先，使用Cutadapt[10]去除测序过程中测序接头、低质量序列及不确定序列，并使用FastQC对所得序列进行质量评估和可信度分析，然后利用Bowtie2[11]和Tophat2[12]软件将测序序列Mapping到捻转血矛线虫的基因组，再使用StringTie[13]做转录本组装，将得到的Clear reads合并，利用perl脚本做转录组从头组装，得到最终转录本。将两个样本Genes的表达量进行FPKM值归一化处理。在差异表达基因(differentially expressed genes，DEGs)的筛选中，使用StringTie和Ballgown来评估所有转录本的表达水平。通过StringTie对两样本基因表达量进行FPKM值归一化，通过Ballgown进行两样本间的差异表达分析。
1.5 差异表达基因注释分析捻转血矛线虫阿苯达唑耐药株给药前后比较，显著差异表达基因定义为P≤0.05且|log2 (fold change)|>1的基因，将获得的差异表达基因分别进行Gene Ontology(GO)功能注释和KEGG富集分析。本研究中，GO功能显著性富集分析首先把所有显著性差异表达基因向GO数据库的各term映射，计算每个term的基因数目，然后应用超几何检验，找出与整个基因组背景比对，在显著性差异表达基因中显著富集的GO条目，通过GO功能显著性富集分析能确定差异表达基因的主要生物学功能。KEGG显著性富集分析以Pathway为单位，应用超几何检验，找出与整个基因组背景比对，筛选出在显著差异表达基因中显著性富集的Pathway，确定差异表达基因参与的生化代谢途径及信号转导途径。
1.6 差异表达基因的荧光定量PCR验证随机选择6个差异表达基因，利用Primer设计引物，同时以β-tubulin为内参基因(相关引物因版面所限未提供)，采用TB GreenTM Premix Ex TaqTM Ⅱ在QuantStudio 7 Flex荧光定量PCR系统进行检测，每个样品设3个重复。经2-△△Ct法数据分析后，利用单因素方差分析数据显著性。
2 结果 2.1 转录组测序对捻转血矛线虫阿苯达唑耐药株给药前和给药后分别进行转录组测序，将得到的数据进行质量控制及组装，结果显示，给药前和给药后分别获得1亿268万以上的原始Reads，且Valid read比例均在99%以上，表明建库质量良好。将两个样本的Genes表达量进行FPKM值归一化处理，如表 1显示，给药前转录本FPKM值＞0有15 354个Genes，FPKM值≥10有2 905个Genes；给药后转录本中FPKM值＞0有14 366个Genes，FPKM值≥10有2 805个Genes；给药前和给药后FPKM值≤10的转录本均占80%左右。
|
表 1 样本有效RNA-Seq序列及差异表达基因数 Table 1 Number of high-quality valid reads and differentially expressed genes of each sample |
选用Ballgown软件对基因表达量做差异表达分析，以显著水平P≤0.05及log2 (fold change)>1或log2 (fold change)＜-1为标准筛选差异表达基因。以log2(fold change倍数变化)为横坐标，-log10(P value)为纵坐标，对差异表达分析中所有的基因绘制火山图(图 1)。将给药前后的差异表达基因进行统计，结果表明，差异表达基因共851个，其中表达基因上调有584个，表达基因下调有267个。表 2列出了差异倍数最大的前10个上调基因和前10个下调基因。此外，与耐药性相关的多数已知基因如ATP酶、糖基转移酶、己糖激酶、谷胱甘肽S转移酶、热休克蛋白等的表达都发生显著变化(表 3)。
|
图 1 给药前和给药后相比差异表达的基因 Fig. 1 Differential expressed gene before and after using albendazole |
|
|
表 2 差异倍数最大的前10个上调基因和前10个下调基因 Table 2 Top 10 up-regulated and down-regulated genes with high differential expression |
|
|
表 3 抗性机制有关的部分重要基因在阿苯达唑给药前后的表达情况 Table 3 Expressions of some important genes related to resistance mechanism before and after using albendazole |
将测序得到的差异表达基因与GO、KEGG数据库进行比对分析，结果显示，经GO功能注释分类后，给药前和给药后共有851个差异表达基因富集到55条GO Term上，其中生物学过程274条，细胞组分117条，分子功能164条。在生物学过程分类中，代谢过程(metabolic process)、翻译(translation)和蛋白水解(proteolysis)等GO Terms富集水平最高；在细胞组分分类中，细胞质(cytoplasm)、细胞膜(membrane)和细胞核(nucleus)等GO Terms富集水平最高；在分子功能分类中，核苷酸结合(nucleotide binding)、ATP结合(ATP binding)和金属离子结合(metal ion binding)等GO Terms富集水平最高。
2.4 差异表达基因KEGG通路分析捻转血矛线虫响应阿苯达唑药物选择压力下显著差异表达基因共851个，对DEGs进行KEGG pathway分析，结果显示，DEGs被富集到的通路有核糖体(ribosome)、吞噬体(phagosome)、细胞凋亡(apoptosis-multiple species)、过氧化物酶体增殖物激活受体信号通路(PPAR signaling pathway)等75条通路上，其中有11条显著富集通路(表 4)。另外与抗性相关的pathway还包括药物代谢——谷胱甘肽代谢(glutathione metabolism)、AMPK信号通路(AMPK signaling pathway)、ABC转运蛋白(ABC transporters)等。
|
图 2 给药前和给药后差异表达基因GO分类注释图 Fig. 2 Gene Ontology classification of the DEGs before and after using albendazole |
|
|
表 4 差异表达基因显著富集的KEGG通路 Table 4 The significantly enriched KEGG pathways in the DEGs |
本试验从转录组测序结果中筛选了差异表达基因，并随机选取6个基因，以β-tubulin为内参基因，利用荧光定量PCR对转录组结果进行验证。结果显示，荧光定量PCR检测的6个基因差异表达倍数与转录组测序表达结果基本趋于一致(图 3)。转录组测序与荧光定量PCR检测的差异倍数存在差异，可能是两种计算方法和试验方法不同所致。
|
**.P < 0.01 图 3 荧光定量PCR验证RNA-Seq试验结果 Fig. 3 Validation of RNA-Seq results by fluorescence quantitative PCR |
近年来，随着高通量测序技术的不断发展，RNA-Seq应运而生，并被广泛应用于生命科学的研究，而寄生虫转录组学测序与分析方兴未艾，为寄生虫耐药性的研究提供了新方法。本试验通过捻转血矛线虫阿苯达唑耐药株给药前后的比较转录组学测序分析，从捻转血矛线虫bBZ和aBZ的RNA-Seq序列信息中鉴定出15 943个差异表达基因，以P≤0.05且|log2 (fold change)|>1为筛选条件，筛选出851个差异显著性表达基因，并对这些差异表达基因在GO和KEGG数据库进行功能预测和分类注释。GO数据库分析表明，Hc_bBZ vs Hc_aBZ注释为特定GO显著差异表达基因数中，细胞组分和分子功能较多，且均在核苷酸结合、ATP结合和金属离子结合，导致捻转血矛线虫对阿苯达唑耐药性产生的基因主要集中在这部分细胞组分和分子功能中，推测其细胞组分和分子功能可能与捻转血矛线虫的耐药性之间存在一定的关系，具体机制尚未清楚，需要在后续试验中深入研究。
通过对差异表达基因的Pathway统计，信号通路主要集中于谷胱甘肽代谢、过氧化物酶体增殖物激活受体(PPAR)信号通路、糖异生/糖酵解和细胞凋亡通路等。谷胱甘肽转移酶(glutathione S-transferase，GST)是生物体内重要的解毒酶，目前在昆虫对杀虫剂耐药性研究较多。Kandil等[14]2017年将GST作为包被抗原，建立绵羊捻转血矛线虫自然感染的ELISA检测方法，结果表明其敏感性和特异性均较好，可用于检测绵羊捻转血矛线虫的感染。本试验发现捻转血矛线虫GST基因表达量在给药前后差异显著，并且通过GO和KEGG分析，注释到谷胱甘肽代谢信号通路，但该基因是否可以用作早期鉴别诊断抗原，有待进一步验证。PPAR是一类配体激活的核转录因子超家族成员，包括PPAR-α、PPAR-β和PPAR-γ 3种表型，其中以PPAR-γ的研究最为深入，现有的研究表明PPAR参与机体脂类代谢、炎症反应和细胞生长等，调节体内糖平衡、能量平衡等多种生物学功能，而捻转血矛线虫在给药前后糖代谢、脂代谢及能量代谢等三羧酸循环相关的基因表达量均差异显著。糖异生/糖酵解代谢途径是生物体提供能量和中间代谢物产生的核心枢纽，是真核生物中普遍存在的保守的基本代谢途径，本试验中能量代谢和糖酵解介导的ATPase在给药后均上调表达。研究表明，伊维菌素能显著影响捻转血矛线虫P糖蛋白高表达，Pgp作为一种能量依赖性药物外排泵，通过ATP供能将细胞内的药物泵出细胞外，降低细胞内药物浓度，从而产生耐药性[15]。本试验通过KEGG数据库分析发现，差异表达基因集中在谷胱甘肽代谢、PPAR信号通路、三羧酸循环等信号通路中，与前人有关耐药性的研究基本一致。
另外，经过筛选分析，从捻转血矛线虫中鉴定出多个热应激蛋白(Hsp)，发现给药后捻转血矛线虫体内的Hsp家族(Hsp20、Hsp70、Hsp90)表达量均上调，并参与多个生物学过程，笔者认为Hsp可能与捻转血矛线虫响应阿苯达唑药物选择压力和对环境的应激有关，在GO数据库中，Hsp家族差异基因均被注释到“对应激的反应”和“防御反应”功能中，其中“对应激的反应”功能注释更加表明了Hsp可能与捻转血矛线虫耐药性有关。Martis等[16]通过对鸡蛔虫氟苯达唑处理前后转录组测序，结果表明，Hsp基因显著上调，此外，在杜氏利什曼原虫、白背飞虱等研究也发现类似结果，这些结果均表明Hsp家族与耐药性有关[17-18]。笔者实验室正在开展差异表达基因的验证及耐药虫株与敏感虫株深入的比较转录组学研究工作，以明确其在捻转血矛线虫对阿苯达唑耐药性中所起的作用，为耐药性早期鉴别诊断方法的建立奠定基础。
4 结论通过转录组测序(RNA-Seq)获得了捻转血矛线虫给药前后的差异表达基因及抗性显著表达上、下调基因，qRT-PCR表达量验证分析表明转录组数据可信度较高。差异表达基因主要富集于ATP结合、代谢过程及蛋白水解等GO条目和谷胱甘肽代谢通路、PPAR信号通路、糖异生/糖酵解通路和细胞凋亡通路，这些基因的相互协同调控可能是捻转血矛线虫对阿苯达唑产生耐药性的重要机制。
| [1] | KAMINSKY R, DUCRAY P, JUNG M, et al. A new class of anthelmintics effective against drug-resistant nematodes[J]. Nature, 2008, 452(7184): 176–180. DOI: 10.1038/nature06722 |
| [2] | MUCHIUT S M, FERNÁNDEZ A S, STEFFAN P E, et al. Anthelmintic resistance:management of parasite refugia for Haemonchus contortus through the replacement of resistant with susceptible populations[J]. Vet Parasitol, 2018, 254: 43–48. DOI: 10.1016/j.vetpar.2018.03.004 |
| [3] | BARRERE V, FALZON L C, SHAKYA K P, et al. Assessment of benzimidazole resistance in Haemonchus contortus in sheep flocks in Ontario, Canada:comparison of detection methods for drug resistance[J]. Vet Parasitol, 2013, 198(1-2): 159–165. DOI: 10.1016/j.vetpar.2013.07.040 |
| [4] | MICKIEWICZ M, CZOPOWICZ M, GÓRSKI P, et al. The first reported case of resistance of gastrointestinal nematodes to benzimidazole anthelmintic in goats in Poland[J]. Ann Parasitol, 2017, 63(4): 317–322. |
| [5] | BARRÈRE V, ALVAREZ L, SUAREZ G, et al. Relationship between increased albendazole systemic exposure and changes in single nucleotide polymorphisms on the β-tubulin isotype 1 encoding gene in Haemonchus contortus[J]. Vet Parasitol, 2012, 186(3-4): 344–349. DOI: 10.1016/j.vetpar.2011.11.068 |
| [6] | ZHANG Z Z, YANG X, AWAIS A A, et al. Development of a tetra-primer ARMS-PCR for detecting the E198A SNP in the isotype-1β-tubulin gene of Haemonchus contortus populations in China[J]. Vet Parasitol, 2018, 252: 127–130. DOI: 10.1016/j.vetpar.2018.01.021 |
| [7] |
薄新文, 李祥瑞. 多重PCR检测捻转血矛线虫苯并咪唑类药物抗性等位基因[J]. 中国农业科学, 2005, 38(4): 826–830.
BO X W, LI X R. Multiplex PCR detection of allele on benzimidazole-resistance or -susceptibity in natural populations of Haemonchus contortus[J]. Scientia Agricultura Sinica, 2005, 38(4): 826–830. DOI: 10.3321/j.issn:0578-1752.2005.04.028 (in Chinese) |
| [8] | LAING R, KIKUCHI T, MARTINELLI A, et al. The genome and transcriptome of Haemonchus contortus, a key model parasite for drug and vaccine discovery[J]. Genome Biol, 2013, 14(8): R88. DOI: 10.1186/gb-2013-14-8-r88 |
| [9] | GASSER R B, SCHWARZ E M, KORHONEN P K, et al. Understanding Haemonchus contortus better through genomics and transcriptomics[J]. Adv Parasitol, 2016, 93: 519–567. DOI: 10.1016/bs.apar.2016.02.015 |
| [10] | MARTIN M. Cutadapt removes adapter sequences from high-throughput sequencing reads[J]. EMBnet J, 2011, 17(1): 10–12. DOI: 10.14806/ej.17.1 |
| [11] | LANGMEAD B, SALZBERG S L. Fast gapped-read alignment with Bowtie 2[J]. Nat Methods, 2012, 9(4): 357–359. DOI: 10.1038/nmeth.1923 |
| [12] | KIM D, PERTEA G, TRAPNELL C, et al. TopHat2:accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions[J]. Genome Biol, 2013, 14(4): R36. DOI: 10.1186/gb-2013-14-4-r36 |
| [13] | PERTEA M, PERTEA G M, ANTONESCU C M, et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads[J]. Nat Biotechnol, 2015, 33(3): 290–295. |
| [14] | KANDIL O M, GAMIL I S, HENDAWY S H M, et al. Efficacy of glutathione-S-transferase purified antigen of the gastro-intestinal nematode Haemonchus contortus in diagnosis of sheep haemonchosis[J]. J Parasit Dis, 2017, 41(4): 968–975. DOI: 10.1007/s12639-017-0920-8 |
| [15] | WILLIAMSON S M, WOLSTENHOLME A J. P-glycoproteins of Haemonchus contortus:development of real-time PCR assays for gene expression studies[J]. J Helminthol, 2012, 86(2): 202–208. DOI: 10.1017/S0022149X11000216 |
| [16] | MARTIS M M, TARBIAT B, TYDÉN E, et al. RNA-Seq de novo assembly and differential transcriptome analysis of the nematode Ascaridia galli in relation to in vivo exposure to flubendazole[J]. PLoS One, 2017, 12(11): e0185182. DOI: 10.1371/journal.pone.0185182 |
| [17] | ZHOU C, YANG H, WANG Z, et al. Comparative transcriptome analysis of Sogatella furcifera (Horváth) exposed to different insecticides[J]. Sci Rep, 2018, 8: 8773. DOI: 10.1038/s41598-018-27062-4 |
| [18] | KAUR P, GARG M, HOMBACH-BARRIGAH A, et al. MAPK1 of Leishmania donovani interacts and phosphorylates HSP70 and HSP90 subunits of foldosome complex[J]. Sci Rep, 2017, 7: 10202. DOI: 10.1038/s41598-017-09725-w |