捻转血矛线虫阿苯达唑耐药虫株给药前后miRNA差异表达谱及网络调控分析

王腾宇 陈昕迪 石雅琴 毛晓伟 闫 旭 苏 娅 温海峰 刘春霞 王文龙*

(1.内蒙古农业大学 兽医学院/农业部动物疾病临床诊疗技术重点实验室，呼和浩特 010018； 2.内蒙古农业大学 生命科学院，呼和浩特 010018)

捻转血矛线虫(

Haemonchus

contortus

)属于毛圆科(Trichostrongylidae)血矛属(

Haemonchus

)，是反刍动物消化道寄生虫中感染性最强的虫种之一，可引起寄生宿主贫血，腹泻等临床症状，严重可导致幼畜及病畜的大批量死亡。目前的主要治疗方法是服用阿苯达唑和伊维菌素等化学药物，由于缺乏科学的指导用药方法，已造成全球多国家和地区呈现多重耐药的局面，给畜牧业带来了巨大的经济损失。为进行科学合理的指导用药，防止耐药性的产生，探究捻转血矛线虫耐药性的产生及调控机制至关重要。

随着基因组学的不断发展，研究发现能够编码蛋白的基因占比很少，其余大部分为非编码RNA。miRNA是非编码RNA中的一类内源性小分子RNA，由18～22个核苷酸构成，广泛存在于动物，植物和微生物中。其主要通过RNA聚合酶Ⅱ，Argonaute-1以及其他转运蛋白发挥调控作用。miRNA是具有高度保守性的单链RNA,可以通过完全或不完全的结合方式与靶向基因的3’非翻译区反应，从而对基因的表达量起到调控作用。在有关昆虫的研究中发现，miRNA主要参与细胞的增殖与分裂，生长发育以及细胞调亡等生物过程。Isik等在秀丽隐杆线虫的研究中发现，hco-miR-34在三期幼虫时期高度表达，说明miRNA参与到机体的应激调节功能。在耐药性的相关研究中，Mortezaei等发现阿苯达唑亚砜对羊的细粒棘球绦虫不同发育阶段具有较强影响，原头蚴时期let-7和miR-61在不同药物浓度下均有显著差异表达。Pérez等研究囊尾蚴对吡喹酮的耐药作用时发现，miRNA全表达谱中大量miRNA表达量出现显著变化。Gillan等发现miR-9551在耐伊维菌素捻转血矛线虫中表达量上调6倍，在抗性圆锥虫中上调16倍。在肺癌等肿瘤细胞中，也存在许多表达量大幅度变化的miRNA直接或间接导致化疗药物的耐药性产生。Marks等通过对捻转血矛线虫微阵列的分析发现hco-miR-228-5p和hco-miR-235-3p在三期幼虫中差异表达量最高，可能在生殖发育等方面发挥重要作用，同时hco-miR-60-3p和hco-miR-236-3p 在四期幼虫中被显著富集，表明其可能参与调节肠道发育。在过去的二十年中，寄生蠕虫在基因组和转录组以及生物信息学等方面取得了重大进展，提高了分子水平上的认识和理解。然而与秀丽隐杆线虫这样的模型动物相比，捻转血矛线虫在转录组学、蛋白质组学和代谢组学等领域的研究仍处于起步阶段，尚未找到能够调节耐药性的关键分子调控组合，仍需要进一步探索。

本研究通过比较捻转血矛线虫阿苯达唑耐药虫株给药前和给药后虫体转录本的差异，探究耐药虫体在受到药物刺激的情况下产生的抗性机制，进一步揭示其耐药机理，为线虫耐药性早期快速检测以及药物作用靶点的开发提供参考。

1 材料与方法

1.1 试验动物与虫株

本研究中捻转血矛线虫阿苯达唑耐药虫株采自内蒙古自治区兴安盟乌兰浩特市科尔沁右翼前旗，选取同群的10只绵羊进行虫卵孵化抑制试验和体内驱虫试验后，通过粪便虫卵检查和感染强度测定等方法确定捻转血矛线虫感染情况并确定其所感染的虫株为耐药虫株。依据阿苯达唑药代动力学特性在第一次给药后25 d剖检2只绵羊采集给药前虫体，随后正常剂量进行阿苯达唑二次给药，在12 h左右血药浓度达到最高时剖检2只绵羊并采集给药后虫体，采集过程中逐条鉴别虫体雌雄并分装于冻存管中，标记好给药前后置于液氮罐中冷冻带回实验室。

1.2 主要试剂

RNAiso Plus(Trizol)，Recombinant DNaseⅠ和PrimeScriptRT reagent Kit with gDNA Eraser反转录试剂盒均购自TaKaRa公司；Mir-X miRNA First-Strand Synthesis Kit(Clontech，美国)，Next Small RNA文库制备试剂盒(New England Biolabs，美国)。

1.3 总RNA提取及质检

分别选取捻转血矛线虫耐药虫株给药前和给药后雌虫和雄虫虫体各1管，每管50条，进行等质量称取和等比例混合，按照RNAiso Plus试剂盒说明书进行RNA提取和纯化，使用Recombinant DNaseⅠ去除RNA样品中基因组DNA，使用NanoDrop2 000和Agilent 2 100检测RNA的纯度以及浓度，高于500 ng/μL的样品且RIN值≥7作为选取标准。

1.4 Small RNA文库构建与测序

将质检合格的Total RNA采用15%变性聚丙烯酰胺凝胶(SDS-PAGE)电泳选取18～30 nt范围的条带切胶，用TRNA连接酶对小片段RNA进行纯化并连接3′和5′接头。运用SuperScriptⅡ反转录酶进行反转录并进行PCR扩增。将PCR产物进行胶回收并纯化140～150 bp的条带，溶于EB solution，完成文库构建。采用Illumina Hiseq2 000平台按照标准化流程完成文库上机测序。

1.5 测序后的数据处理与miRNA靶基因预测分析

运用FastQC软件对下机的原始数据(Raw data)进行质量评估，过滤掉数据中的低质量reads，不带有3′接头以及含有5′接头的reads和含有polyA的reads，最后筛选出大于18 nt且小于30 nt的reads，将得到高质量的clean reads进行下一步分析。对测得的Small RNA进行Rfam数据库注释，并将clean reads通过Bowtie软件比对捻转血矛线虫参考基因组(GenBank：GCA_000 469 685.2)，确定tags在参考基因组的具体位置。

使用软件Repeat Masker version open找出repeat associate的Small RNA序列，对已存在的miRNA和已知的miRNA进行鉴定，确定miRBase数据库与参考基因组的比对率。对得到的样品中已存在miRNA的含量进行统计，比对miRNA前体序列，过滤去除鉴定为来自mRNA降解片段，repeat区域，以及rRNA、scRNA、snRNA和tRNA等其他小RNA的tag序列，挑选出能够比对上参考基因组的tags。

使用MIREAP预测miRNA特殊的二级结构，并根据miRanda的算法将未能匹配上的序列进行Dicer酶切位点及能量等特性分析，预测出可能存在的新miRNA。最终利用miRNA种子区匹配保守性位点，使用RNAhybrid、Miranda和TargetScan 3个生物信息学系统进行靶基因预测和评分估计，根据3种不同算法并结合二聚体热力学稳定性等特征，将靶基因预测结果取交集作为miRNA靶基因预测的结果。

1.6 差异表达miRNA生物信息学分析

对耐药虫株给药前和给药后虫体的miRNA进行差异表达分析，将edgeR软件dispersion 设置为0.01，其他设为默认参数，表达量变化2倍以上并且FDR值<0.05作为筛选的标准阈值。为进一步识别给药前后显著差异表达的miRNA，通过计算miRNA的TPM值对其表达量进行归一化处理，公式为TPM=T*10/

N

(T表示miRNA的tags，

N

表示总miRNA的tags)。将DE miRNA所对应的靶基因应用KOBAS软件向GO数据库的GO terms进行映射，同时进行靶基因的KEGG Pathway显著性富集分析，经过FDR校正后，选择

Q

-value≤0.05的Pathway定义为候选基因中显著富集的Pathway。对注释富集结果以分类二级柱状图和圈图的形式展现不同差异目的基因集合在GO terms和KEGG Pathway的分布情况，重点关注可能参与抗性发生的代谢调控通路。

1.7 DE miRNA靶向DE mRNA关联分析

由于miRNA的5’端2～8个碱基可与其靶向的mRNA的3’-非翻译区完全匹配互补，表达调控模式呈现负相关，为深度解析耐药虫株在药物影响下DE miRNA参与调控的基因，将筛选出的DE miRNA与其趋势相反的mRNA建立调控关系互作网络，使用Cytoscape软件对DE miRNA的靶向网络进行可视化作图。

1.8 转录组DE miRNA Stem-loop qRT-PCR

为保证测序数据的可靠性，随机挑选5个DE miRNA采用Stem-loop qRT-PCR验证表达量调控趋势，每组样品设置3个生物学重复以及技术重复，根据5个DE miRNA的成熟体序列利用DNAMAN软件设计特异性的茎环引物及qPCR引物，委托生工生物工程(上海)股份有限公司合成，引物序列见表1。将Trizol法提取的捻转血矛线虫总RNA分别利用特殊设计的Stem-loop引物进行反转录，以每个miRNA相对应的cDNA为模板通过ABI QuantStudio7 Flex system进行qPCR测定，反应条件为94 ℃预变性30 s，95 ℃变性5 s，60 ℃延伸35 s，共进行40个循环。选取U6 snRNA作为内参基因，最后采用2计算相对表达量。

表1 DE miRNA序列及Stem-loop qRT-PCR引物
Table 1 DE miRNA sequences and Stem-loop RT-qPCR primers

miRNA引物序列(5′-3′)Primer sequence (5′-3′)hco-miR-5960RT primer: GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACATAACCF: AACTTGATGTGGTCCGGAGTChco-miR-5885cRT primer: GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTAGCGAF: TCCAGCTGAGATCACGCGTAhco-miR-259RT primer: GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTGCAACF: TCGCGTAATCTCATCCGAATCTGhco-miR-228RT primer: GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCCGTGAF: TCCGAAATGGCACTGCATGAAnovel-m0498-5pRT primer: GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCTCACAF: AACCACTCGTTGTAGGACCCU6F: GAGAAGATTAGCATGGCCCCTR: AAAATTTGGAACGCTTCACGAA通用下游引物Unified reverse primerGTCGTATCCAGTGCAGGGT

2 结果与分析

2.1 RNA质检分析

经Nanodrop确定指导区间后，Agilent检测给药前耐药虫株和给药后耐药虫株总RNA浓度分别为498和510 ng/μL，在18 S和28 S位置均出现波峰，RNA完整性好，质量评级为A类。

2.2 Small RNA数据评估及长度分布统计

从捻转血矛线虫给药前和给药后2个对照组中分别获得14 925 456和15 304 062条Total reads，过滤掉低质量序列和去接头后得到13 888 536和12 928 074条高质量可比对序列，占比94.623 8%和86.036 9%，测序质量良好，数据准确可靠，可以进行后续分析和试验(表2)。对miRNA长度分布进行统计，大多数序列的长度主要集中在20～22 nt，在22 nt处频率最高出现波峰(图1)。将给药前与给药后2个样本文库经过滤后的clean reads分别比对到GenBank和Rfam数据库，最大程度上找到并除去样本中可能的rRNA、scRNA、snoRNA、snRNA和tRNA。比对GenBank结果显示，rRNA分别占8.02%(1 113 201)和13.97%(1 806 675)，tRNA分别占0.01%(1 520)和0.10%(12 984)，比对到Rfam数据库的结果显示，rRNA分别占1.12%(155 494)和1.82%(234 689)，snRNA占0.01%(1 627)和0.03%(3 233)，tRNA分别占0.21%(29 165)和0.30%(38 336)，snoRNA分别比对到148条和202条tags。

2.3 miRNA分类鉴定

通过比对miRBase数据库对不同序列的丰度进行统计，给药前测得的序列中，miRNA占83.59%；给药后的序列中，miRNA占62.24%(图2)。其中耐药虫株给药前与给药后分别鉴定到已公布的446个和553个已存在的捻转血矛线虫miRNA；分别鉴定到133个和125个没有比对到捻转血矛线虫miRNA的tags序列，但在其他物种中已知的miRNA；同时鉴定到876个和659个novel miRNA。此外，对耐药虫株给药前和给药后已知的mRNA进行了碱基偏向性分析， Dicer酶在识别和切割前体miRNA时，5′端首位碱基对U碱基具有很强的偏向性，在长度集中的18～25 nt的miRNA中，两样本的首位碱基均主要偏向于U，其次偏向于A和C(图3)。

表2 数据预处理与去端头情况统计
Table 2 Data filtering and de-joining statistics

样本Sample总读值Total reads3′端缺失Readswith null3′ adapter5′端污染的读值5′ adaptercontaminants小于18 nt的读值Reads shorterthan 18 nt含有poly(A)的序列Reads withpolyA低频数读值(<2)Readswith Lowfrequency干净读值Cleanreads给药前虫株Pre-drug strain14 925 456100 5317 114327 308186338 09913 888 536给药后虫株Post-drug strain15 304 062171 56724 9501 243 235121635 70812 928 074

图1 给药前后Small RNA长度分布Fig.1 Small RNA length distribution before and after administration

图例括号中的数字表示该条目tags数占总tags的比值。 Numbers in brackets in legend represent the ratio of the number of tags in the entry to the total tags.图2 给药前(a)和给药后(b)不同分类tag丰度统计Fig.2 Abundance statistics of different classification tags before (a) and after (b) drug administration

图3 给药前(a)和给药后(b)特异表达的已知miRNA碱基分布Fig.3 Distribution of known miRNA bases specifically expressed before (a) and after (b) administration of the drug

2.4 DE miRNA的差异表达分析

通过DEseq2软件基于负二项分布的分析，根据read counts计算每个miRNA的TPM表达量，成功获得捻转血矛线虫耐药虫株给药前虫体和给药后虫体的miRNA序列及其表达谱。基于对miRNA表达谱预测结果的比较分析共筛选到125个上调表达的miRNA和63个下调表达的miRNA(图4)。其中miR-10512-y、miR-7630-y、miR-3605-y以及miR-3801-x等miRNA显著上调表达，log(Fold change)分别达到了11.04、9.92、9.03和8.34；novel-m0934-5p，hco-miR-5942、novel-m0110-5p和novel-m0190-5p等miRNA显著下调，log(Fold change)分别达到-6.91、-6.43、-6.42和-6.43。表3为表达量差异表达较高的10个miRNA。

2.5 miRNA靶向mRNA预测及功能富集分析

2

.

5

.

1

miRNA靶基因预测

图4 给药前后组间差异miRNAs火山图Fig.4 Volcano plot of differential miRNAs between groups before and after drug administration

通过对miRNA靶基因的预测和统计分析，由表4可知，给药前耐药虫株1 454个miRNA预测到的靶基因数量为5 919个，对应133 723个靶基因结合位点；给药后耐药虫株1 313个miRNA预测到的靶基因数量为5 863个，对应113 058个靶基因结合位点。

2

.

5

.

2

miRNA靶基因GO和KEGG富集分析

给药前和给药后虫体差异表达的miRNA靶基因GO功能注释显示，miRNA靶基因显著富集到了46条GO条目，其中生物过程24条，细胞组分13条，分子功能9条。由图5可知，GO:0 009 987(细胞过程)、GO:0008 152(代谢过程)和GO:0065 007(生物调节)等生物过程；GO:0044 464(细胞部分)、GO:0043 226(细胞器)和GO:0016 020(膜结构)等细胞组分；GO:0003 824(催化活性)、GO:0005 488(结合)和GO:0005 215(转运活动)等分子功能富集水平最高。

同时对2个比较组进行KEGG富集分析，结果显示差异表达的miRNA靶基因主要富集在碳水化合物代谢通路中的糖酵解/糖异生1013个Genes(ko00010)，氨基酸代谢通路中的半胱氨酸和蛋氨酸代谢238个Genes(ko00270)，氨基酸代谢通路中的缬氨酸、亮氨酸和异亮氨酸降解172个Genes(ko00280)，运输和代谢分解通路中的吞噬体155个Genes(ko04145)等信号通路(图6)。在对通路中富集到的基因数进行统计发现，在基因数大于50的富集通路中，ABC转运蛋白信号通路涉及差异的miRNAs靶基因有117个(ko02010)，HIF-1信号通路涉及差异的miRNAs靶基因有96个(ko04066)，cGMP-PKG信号通路涉及差异的miRNAs靶基因有77个(ko04022)以及细胞色素P450药物代谢通路涉及差异的miRNAs靶基因有66个(ko00982)。

表3 给药前后差异表达较高的10个miRNAs
Table 3 The 10 up- and down-regulated miRNAs with higher differential expression before and after drug administration

miRNA名称miRNA name折叠倍数log2 (Fold change)FDR值FDR value序列Sequencehco-miR-59601.778 284 6520GTGGUCCGGAGUCGGAGGGUUAUnovel-m1107-5p3.850 483 4513.014 939 229UAGGCUUCGGCUUAGGCUUCmiR-9277-y3.753 792 1832.968 520 151CGGUUCGAAUCCUGCGGAmiR-11987-x2.298 296 3735.301 084 290GAGGAAACUCUGGUGGAAGnovel-m0228-5p2.064 840 1671.752 001 699GUGGUCCGGAGCCAGUGGGUGAGUnovel-m0317-3p-2.29 716 8875.620 018 481ACAUGCCAUGCAUCCGGCCGUAhco-miR-5931-1.693 676 3890.000 727 566UGCAGGUCUUGAUUGGAAGGCUmiR-5895-x-1.349 773 6269.674 127 840UACCGUAGCAUCUGUAUGUCnovel-m0498-5p-1.056 174 5763.860 652 208CGUUGUAGGACCCCCUGUGAGhco-miR-40e-1.009 179 9930.000 011 123ACUGCAGGACCGCUCAGGUGAAG

表4 给药前和给药后miRNA靶基因位点统计
Table 4 Predictive statistics of miRNA target genes locibefore and after drug administration

样本名称Simple namemiRNA数量miRNA number靶基因数量Target gene number靶基因位点数量Target site number总数量Total number2 76711 782246 781给药前虫株Pre-drug strain1 4545 919133 723给药后虫株Post-drug strain1 3135 863113 058

图5 给药前后DE miRNA的靶基因GO数据库注释分类图Fig.5 GO database annotated classification chart of target genes of DE miRNA before and after drug administration

2.6 DE miRNA互作调控网络的构建

结果显示144个DE miRNA共结合到143个基因，共涉及373个miRNA-mRNA调控关系组合，其中92个上调的miRNA共靶向到178对miRNA-mRNA组合，52个下调miRNA靶向到了195对miRNA-mRNA组合。由图7的可视化分析可以看出同一个miRNA可以靶向结合多个mRNA，同一个mRNA也可以靶向多个miRNA进行调控。

图6 给药前后DE miRNA的靶基因KEGG数据库注释圈图Fig.6 Annotated circle of KEGG database of target genes of DE miRNA before and after drug administration

表5 给药前后与耐药相关的显著差异表达miRNAs
Table 5 Significantly differentially expressed miRNAs associated with before drug resistance and after drug administration

功能描述Function description小RNAmiRNA类型Type基因名称Gene name类型Type折叠倍数log2 FCFDR值FDR value含有热休克蛋白 Hsp20 结构域和转酮醇酶结构域的蛋白质Heat shock protein Hsp20 domain and Transketolase domain containing proteinmiR-234-y下调HSP20PDHB-1上调1.746 2431.513 9313.8×10-57.58×10-35含有蛋白质的丝和中间丝结构域Filament and intermediate filament domain containing proteinhco-miR-5960上调IFA-4下调7.32 1931.02×10-3兰尼碱受体Ryanodine receptormiR-10023-x上调RyR下调8.02 2377.34×10-19UTP-1-磷酸葡萄糖尿苷转移酶和含有蛋白激酶 C 结构域的蛋白质UTP-glucose-1-phosphate uridylyltrans-feraseand Protein kinase Cdomain containing proteinmiR-3605-x上调UGP2PKC-1IFA-4下调2.696 422.62 7097.321 932.06×10-148.43×10-91.01×10-3GNS1/SUR4家族蛋白和激酶结合蛋白包含CGI-121结构域的蛋白GNS1/SUR4 family protein and Kinase binding protein CGI-121 domain containing protein novel-m0498-5p下调ELO-5TPRKB上调1.024 2661.890 0821.01×10-481.21×10-25

图7 给药前后与耐药相关差异表达miRNA-mRNA调控网络Fig.7 Differentially expressed miRNA-mRNA regulatory network associated with drug resistance before and after drug administration

2.7 DE miRNA Stem-loop RT-qPCR验证sRNA-seq结果

如图8所示，5个DE miRNA在比较组中的变化情况均通过RT-qPCR结果在组内真实表达出来，上(下)调表达趋势与Small RNA测序结果一致，进一步证明sRNA-seq的准确性和可靠性。

图8 miRNA荧光定量PCR验证测序结果Fig.8 The result of qPCR compared with miRNA-Seq

3 讨 论

随着转录组学和测序技术的快速发展，越来越多的转录本数据和Small RNA数据被用来分析宿主引起线虫感染时基因发生的差异表达。通过对高通量测序的下机数据进行过滤和比对可以得到大量具有研究意义的基因序列，从而发现许多微卫星分子标记。miRNA作为一种调节某些基因转录后时空表达的调控分子，也是重要的表型调控位点，而耐药性的产生可能会影响某些特定miRNA的显著差异表达。大量转录因子不再通过翻译后进行调节，而是在被抑制转录后通过某些途径诱导miRNA辅助细胞限制体内和体外的表型产生，为疾病的治疗和研究提供了新的研究方向。Lucienne等发现犬心丝虫的三期和四期幼虫会分泌大量候选miRNA，由此可以作为新型的诊断标记。Seifollah等通过对细粒棘球绦虫进行阿苯达唑给药发现，2个挑选出的miRNA在不同药物浓度和虫体发育时期均出现了显著差异表达，进一步说明miRNA表达会受到药物浓度的影响且十分灵敏。同时miRNA广泛参与机体早期的胚胎发育，细胞代谢与凋亡等生理过程。由此猜测，miRNA可能通过抑制靶基因的表达从而对捻转血矛线虫的耐药性产生影响。

根据表达量分析可以看出，在药物刺激的作用下通过上调miRNA从而抑制靶基因表达的调控模式占据了主导地位。Pentimone等在研究厚垣孢普可尼亚菌接种番茄根部不同天数miRNA的差异变化时发现，参与宿主对真菌反应的Small RNA在24 nt的部分丰度出现了明显减少。根据本研究的Small RNA长度分布统计显示，miRNA均在21 nt和22 nt分布频率较大且给药前分布数量较多，但小于21 nt和大于22 nt的miRNA数量在给药前则明显少于给药后。耐药虫株在药物作用下出现这种显著性差异，由此推断miRNA在结构特征等方面发生了改变，预示着耐药虫株在药物作用的胁迫下可能通过改变miRNA结构影响耐药性。研究发现在生物体受到侵染的情况下，部分基因会受到不同程度的抑制或激活，且大量miRNA通过靶向各类代谢及免疫途径对受到的胁迫作用产生应答。Kaiser等在cGMP信号通路中发现

PKA

和

PKG

同源基因通过调节环核苷酸水平的环化酶和磷酸二酯酶对疟原虫的耐药性产生调节作用。在本研究中显著差异表达的靶基因也关联到cGMP通路，该通路中

PKG

同源基因

PRKCE

被显著富集。

PRKCE

作为一种蛋白激酶在胆囊癌的研究中与多个miRNA呈负调节，进而对吉西他滨的耐药性产生影响。由此可以猜测

PRKCE

基因是药物调节通路cGMP中的重要调节因子。通过与前期敏感虫株miRNA表达谱的比较发现，

RAC1

和

FZD4

作为miRNA的靶基因在耐药虫株中也被大量富集到Wnt和cAMP等重要的信号传导通路，说明这2个基因在耐药性的调节中发挥巨大作用。通过对捻转血矛线虫给药前后差异表达的miRNA-mRNA网络进行关联分析发现，大量miRNA的上调导致蛋白激酶，锌指蛋白和微管蛋白等结构域蛋白相关基因被抑制表达，另有含细胞色素P450结构域蛋白，P-糖蛋白和免疫球蛋白等相关基因在miRNA被抑制的情况下被其激活并发挥作用。Tan等在罗氏沼虾先天免疫防御细菌和病毒的研究中发现，miR-993是关键的一个代谢调节小RNA，对甲壳纲动物的免疫防御功能发挥重要作用。本研究结果中miR-993的表达水平在药物作用的影响下持续下调，其靶基因

HSP20

和

CG9801

被激活。其中

HSP20

已被报道小菜蛾抗杀虫剂品系较敏感品系出现了较弱的上调。热应激蛋白(HSP)在GO和KEGG富集分析中也被广泛注释到了应激和防御调节，与本研究结果一致。

HSP20

还靶向到了miR-234、miR-137、miR-155和miR-4 703等多个miRNA，在耐药性的调控机制网络中未见有具体功能的相关报道，可将其作为候选调控因子进一步研究。Hoy等发现3个miRNA(bantam、miR-277和mir-3479-3p)可作为感染曼氏血吸虫的直接标志或是可能引起病理性感染的间接标志。在本研究结果中显示miR-227在给药后表达呈现差异上调，抑制了靶基因

UNC-6

的表达，而

UNC-6

在秀丽隐杆线虫中被作为诱导保守乙酰胆碱神经递质的关键因子，帮助调节神经系统发育，乙酰胆碱类激酶和受体是寄生虫神经系统调节的重要调控元件，由此推测

UNC-6

可通过miR-227的上调表达抑制其对乙酰胆碱递质信号的传递，从而在药物选择和外界环境的压力下对信号传导进行干预。本研究在对KEGG通路的统计分析中发现，大量的miRNA和靶基因参与代谢、转化、轴突传导和信号传递等生物学过程，其中有24.32%的基因被显著富集到了代谢途径，10.94%注释到了细胞自噬和降解以及50.45%注释到了各类信号传导通路，因此在药物的影响下miRNA及其靶基因很有可能通过鞘脂信号通路，PI3K-Akt 信号通路和FoxO信号通路等一些重要的信号传导通路参与耐药性的调控。有关捻转血矛线虫阿苯达唑耐药相关的信号通路和调控基因仍未见明确报道，依据获得的捻转血矛线虫耐药虫株给药前后差异表达谱，后续可对筛选出的miRNA及其靶基因与其他非编码RNA进一步关联，并对其所处的代谢通路进行生物学功能验证，从而为阐述在非编码RNA水平研究捻转血矛线虫耐药性提供可能的方向。

4 结 论

通过对捻转血矛线虫给药前后的耐药虫株进行Small RNA测序研究，预测到大量高度同源的miRNA和新miRNA，构建miRNA及其靶基因调控网络发现，这些非编码RNA所构成的调控模块广泛参与代谢和信号传导过程，为进一步揭示miRNA在捻转血矛线虫耐药性调控中的关键作用提供了新的视角，并为耐药性早期分子诊断方法的建立提供了数据资源。